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INTRODUCTION 

This chapter reviews the mechanics of shear slippage and 
rupture in rock masses. The physical problem often arises because 
of the presence of a plane or thin zone of weakness - sometimes a 
joint or a fault - in a body under applied load. Discussion of 
the formation of such a plane of weakness is not within the scope 
of this article. Rather we shall focus on the slip response under 
increasing gross shear load transmitted across a pre-existing weak 
plane. The slip response may be stable initially, then becomes 
unstable, leading to a dynamic shear rupture at higher applied 
loads. When the slip distribution is very non-uniform with slip 
confined to a finite segment of the plane of weakness, a Griffith 
crack type failure may result. In contrast, the whole plane may 
cut through the body and slip occurs more or less uniformly. The 
former type of crack failure is usually described by fracture 
mechanics, and the latter by classical strength theory. A more 
general concept, known as the slip -weakening model (and 
tension- softening model in the corresponding tensile failure mode) 
recognizes the elastic brittle crack description and the simple 
strength description as limiting cases. The connection between 
these concepts forms a major focus of this review article. 

While the applications and examples chosen to illustrate the 
theories presented are heavily slanted towards earthquake faults, 
the concepts are equally valid and applicable in many other fields 
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of research, such as in geotechnical engineering where failure of 
jointed rock masses or failure in overconsolidated clay slopes are 
of major concern. The relevance of the mechanics of slip rupture 
to understanding the physics of the earthquake source process is 
obvious. For an appreciation of the physical context of slip 
rupture representation of earthquakes , the reader is referred to 
the text (especially Chapter 7) by Kasahara (1981) on the relation 
between tectonics and earthquakes. 

Elastic dislocation theory has proved to be very useful in 
describing displacement discontinuities (slip) in otherwise 
continuous bodies . A convenient way of thinking about a 
dislocation is to imagine a planar cut in a body. The two faces 
of this cut are slid uniformly and then are glued back together. 
This operation introduces dislocation stress and displacement 
fields in the body. A non-uniform slip distribution can be 
constructed by putting a number of such dislocations at logistic 
locations, and the resulting stress and displacement fields are 
obtained by superimposing the stress and displacement fields of 
each discrete dislocation. This construction is sometimes 
referred to as a continuously distributed or smeared dislocation. 
Geophysicists have taken advantage of such a procedure to estimate 
the seismic slip distributions on earthquake faults by inverting 
the measured displacement field on the ground surface (see, e.g. 
Chinnery, 1961, 1970; Savage and Hastie, 1966; Walsh, 1969 and 
Rybicki, 1986). More recently, post-seismic leveling measurements 
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have also been used to estimate rheological properties based on 
dislocation displacement fields derived from layered 
elastic/viscoelastic half -space models (see, e.g. Thatcher et al, 
1980) . Such kinematic models occupy an important place in using 
measurable quantities of surface deformation to deduce (directly) 
unmeasurable quantities of fault slip and material properties of 
the earth. However, a shortcoming of kinematic modelling is that 
it provides limited insight into the physical processes leading to 
an instability. Such insights are important, for example, to 
forecast an earthquake rupture, since it would provide a rational 
basis to interpret seismic (such as foreshocks) and aseismic (such 
as surface displacement rate and strain rate changes) precursory 
signals. The study of non-kinematic models of slip rupture, of 
which elastic brittle crack theory is a special case, is another 
major theme in this article. 

To reach beyond kinematic modelling, it is necessary to 
prescribe some kind of constitutive law which relates fault slip 
to stress. In order for a seismic event to radiate energy from 
the source, the stress is expected to drop with increasing slip. 
This implies that natural faults can be described by a 
si ip -weakening model. In the laboratory slip-weakening behavior is 
directly observable in a variety of specimens, including 
overconsolidated clay samples, and in intact, sawcut and jointed 
rock specimens . Extrapolations of material parameters obtained 
from the laboratory to the real Earth has always been a difficult 
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task, and those related to the slip-weakening model are no 
exception. Even so, some recent non-kinematic models (e.g. 
Stuart, et al 1985) have generated results consistent with 
available geodetic data covering a period of time. Moreover, such 
models project into the future how slip continues to develop 
until instability. 

This article is organized in three sections. The first deals 
with a summary of essential ideas in fracture mechanics, 
emphasizing the interpretation and relation among the fracture 
parameters K, G and J in shear cracks. This section is concise 
because of the widely available literature on this subject and 
several recent review articles on their applications to 
geophysical problems (Rudnicki, 1980; Rice, 1980; Dmowska and 
Rice, 1986). The second section describes the slip -weakening 
model.. The physical interpretation of the si ip -weakening model 
and connections to G and J are emphasized. The model is used to 
illustrate the loss of stability of a simple slip system. This 
section also summarizes fracture resistance properties deduced 
from laboratory tests and from observations of earthquake 
faulting. The third and last section deals with the general 
formulation of the problem of non-uniform slip distribution in a 
continuum. There are two focuses in the section: the structure 
of the stress transmission Green's function which incorporates 
information about the rheology (elasticity, viscoelasticity, and 
poro- elasticity) and geometry of the continuum containing the slip 
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plane; and the formulation of non-kinematic problems. Several 
examples from recent geophysical literature are discussed. 

The coverage in this article is necessarily incomplete, even 
as far as the mechanics of shear rupture applied to earthquake 
faults is concerned. For example, the dynamic 'rupture process 
important to seismology is not included. Excellent review of this 
subject could be found in Dmowska and Rice (1986), and in the text 
by Aki and Richards (1980) on quantitative seismology. With few 
exceptions, most of our discussions will focus on 2-D problems. 
This should not be construed as an indication that shear rupture 
problems are inherently 2-D, although in many circumstances they 
could be approximated as such. In section 4, we do describe a 
line-spring procedure which reduces a 3-D problem into a 2-D 
one. Within limitations, such a technique appears to be quite 
powerful and provides a computationally economical alternative to 
solving full-scale 3-D fault problems. 

Apart from the references mentioned above, the reader will 
find the following publications of particular interest: Journal 

of Geophysical Research special issue on Fault Mechanics and its 
Relation to Earthquake Prediction (Vol 84, May, 1979), Pure and 
Applied Geophysics special issue on Earthquake Prediction (Vol 
122, No. 6, 1984/5), and an American Geophysical Union publication 
Earthquake Source Mechanics published for the 5th Maurice Ewing 
Symposium (1986) . 
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II SHEAR FRACTURE MECHANICS 

In contrast to applications to technological materials which 
deal mostly with tensile fracture, application of fracture 
mechanics to the earth's crust involves mainly shear cracks. This 
is due to a number of reasons, notably the presence of lithostatic 
pressure which reduces the tendency of tensile cracking on a large 
scale. Repeated ruptures on the same plane also prevents fracture 
propagation from deviating markedly from the pre-existing weakened 
fault plane. Pre-existing fault zones must act as guides to shear 
rupture since without such a weakened plane, shear cracks in 
brittle laboratory specimens often tend to develop kinking or 
wing cracks from the ends of a shear sliding plane. While 
tensile cracking is not usually seen on a scale of tens or 
hundreds of km, it can still form on a more local scale, sometimes 
in the form of en-echelon tensile cracks which are joined together 
by the through -running main shear rupture. As an example, soil 
cracking at an angle to the main rupture was observed in the 1966 
Parkfield earthquake (Allen and Smith, 1966). Recent work on 

compression failure of brittle material by Nemat-Nasser and Horii 
(1982) and Ashby and Hallam (1986) may shed some light on this 
phenomenon. However, our focus in this article will be on a much 
larger scale than that of the en-echelon tensile cracks. 

In general, there is not much difference between the analysis 
of shear and tensile cracks. Two points which play an important 
role in understanding earthquake ruptures, however, should be 
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noted. The first is that unlike tensile cracks, shear crack faces 
are not stress-free even if the surfaces have been well slid. 
However, the residual frictional resistance may be regarded as a 
reference stress level, as will be explained further in section 
2.1. In any case it is useful to keep in mind that the frictional 
work acts as an additional energy sink, reducing the amount of 
energy available to drive the crack tip. The other characteristic 
in shear cracks is that the fracture energy release rate, at least 
for laboratory rocks , is often two orders of magnitude higher than 
that for tensile cracks (Wong, 1982a) . This is presumably related 
to the difference in the physical break-down processes at the 
crack tip of a shear and a tensile crack. The shear breakdown 
process may involve the extension and linking of smaller scale 
en-echelon tensile cracks as described earlier. Laboratory 
measurements of the critical energy release rates in rocks must 
therefore be made in the fracture mode appropriate to the field 
conditions, although we shall see in section 3 that even this does 
not fully account for the discrepancy between measured magnitudes 
in the laboratory and magnitudes estimated from field 
observations. Other than the differences pointed out above, the 
analysis of shear and tensile cracks are rather similar. In the 
rest of this article, the discussions will focus mainly on shear 
cracks . 

Two well-known approaches have been used in the study of 
elastic brittle cracks. The first is based on the work of Irwin 
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(1960) who characterized the intensity of the crack tip stress 
field by a stress intensity factor (i = I, II, III denoting the 
3 modes of deformation) . It is hypothesized that when reaches 
a certain value, K. , crack extension occurs. This fracture 

1C 

criterion = K^ c represents a balance of the crack driving 

stress intensity with a critical stress intensity (fracture 

toughness) that the material can sustain. The second approach is 

an extension of the idea of Griffith (1920) by characterizing 

fracture as a balance between available energy G to drive the 

crack and the energy absorbed by the inelastic -breakdown processes 

of the material at the crack tip G . This fracture criterion 

c 

G = G c may be shown to .be consistent with = K^ c quasistatic 
crack propagation analysis in a linear elastic body. These 
single -parameter fracture characterizations are analogous to the 
classical strength concept which relates the shear stress a to the 
shear strength a at failure of a specimen deforming uniformly. 
There is, however, a major difference between fracture mechanics 
and the strength concept: The strength concept in general cannot 

characterize objectively the crack driving force and therefore 
fails to predict the load level a structure with flaws can carry. 
In section 3, it will be seen that the strength concept and 
brittle elastic crack mechanics may be regarded as two opposite 
limiting conditions of a more general slip -weakening model. 

In the following, we shall summarize the essentials of 


I 
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elastic brittle crack mechanics with particular focus on the 
relationships between the various fracture parameters K, G and the 
path- independent J- integral. The J- integral extends the realm of 
linear elastic fracture mechanics (LEFM) to situations where the 
small scale yielding condition (to be explained) of LEFM is 
violated. 

2 . 1 Elastic Brittle Crack Mechanics 

2.1.1 Asymptotic Crack Tip Stress Field: Near a sharp crack 

such as that shown in Fig. 1, the stress field based on a linear 
elastic analysis may be expressed in the asymptotic form (see, 
e . g. , Rice , 1968a) . 

For mode II deformation, 

ct 12 “ K II (27rr)* 1/2 cos (0/2) [1 - sin(0/2) sin(30/2)] +a f + o(r 1/2 ) 

°22 ” K II ( 2 *r)' 1/2 sin(0/2) cos(30/2) + + o(r 1/2 ) (la) 

CT 11 “ K II (2^r)' 1/2 sin (0/2) [2 + cos(0/2) cos(30/2)] + o(l) 

and the near tip crack face shear slip (for plane strain) is given 
by 


*i - < - "i - -^Hr 1 K n < 8r/ ” > 


. .1/2 A . 3/2 
r/n) ' + o(r ' 


(r ' ) (lb) 


where u^ and u^ are the displacements in the x^ direction for the 
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upper and lower crack faces , and n and v are the shear modulus and 
Poisson's ratio respectively. For plane stress deformation, the 
factor (l-i/ ) should be replaced by l/(l+i/) . 

For mode III deformation, 

a 13 ” - K III (2^ r )' 1/2 sin(0/2) + o(l) 

° 22, “ K m (2*r)' 1/2 cos ( 6 / 2 ) + a f + o(r 1/2 ) (2a) 

and the near tip crack face shear slip is given by 

S 3 = u+ - u^ - (K iiz /m) (8rA) 1/2 + o(r 3/2 ) ( 2 b) 

The first terms in the stress expressions are 1 /,/r singular, 
dominating the behavior of the stress field near the crack tip, 
with singularity strength given by the stress intensity factor K. 
(For brevity, K will stand for K.^ and and S will stand for 

5 ^ or 63, and should be clear from the context of discussion). 

f 

The constant terms a and are retained to show explicitly the 
possibility of non-zero tractions on the crack faces and these 
terms represent limits of shear and normal tractions as the tip is 
approached from within the crack. Clearly no real material can 
withstand infinite stresses so that at the crack tip, the material 
behaves inelastically (shown schematically as the lightly shaded 
region at crack tip, Fig. 1 ). The continuously rising <7_ as 
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r -+ 0 is therefore an artifact of the assumption of elastic 
behavior in the analysis from which (1) and (2) are derived. At 
larger distances the terms left out .in the asymptotic expansion 
become significant and should not be ignored. The stress field 
(la) and (2a) are therefore valid in an annular region (shown as 
shaded in Fig. 1) surrounding the crack tip. This is often 
described as the K- dominated region. 

Equation (1) and (2) completely describes the spatial 
distribution of the near- tip stress field and the crack face 
displacement. This form is independent of the geometry of the 
body containing the crack, and does not depend on the particular 
manner in which the body is loaded. This information is contained 
in K, which is indeed the utility of (1) : For any geometry and 
loading where K is known, the near tip stress field and crack face 
displacement are completely specified. Since K defines the 
strength of the stress singularity at the crack tip, it is the 
parameter which is used in Irwin's fracture criterion mentioned 
earlier . 

Many solutions have been obtained for elastostatic crack 
problems. For a summary of solution methods, see Parker (1981). 
Solutions for the stress intensity factor K are tabulated in Tada 
et al (1973) and Paris & Sih (1965). It is useful to recognize 
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(by dimensional considerations) that K must have the form 

K - aj L F(geometry, loading) (3) 

where a is a generalized loading stress and L is a characteristic 
length in the geometry of the body (often the crack length or half 
crack length). The non-dimensional function F depends on the 
details of geometry and loading. For example, for a center crack 
with half crack length i in a plate loaded remotely by o° and with 
uniform shear resistance a ^ (Fig. 2), 

= ( a°-o ^) Ji Jn or 

•= A a Jk$. where A a = (o° -a^) 

showing that F is a constant (Jir) in this simple case. The second 

form (4b) suggests the usual interpretation of stress drops A a in 

seismological literature. In that case, a° is the shear stress 

f 

prior to seismic rupture and a is the residual friction on the 

ruptured fault segment. The in- situ absolute values of the stress 

states a° and o ^ are not readily determinable. However, it is 

usually the stress (or strain) change that is of interest. Hence 
f 

a may be regarded as a reference stress state. 

Another example is a semi - inf inite crack in an infinite body 
whose crack faces are loaded by line forces P at a distance b from 
the crack tip (in force/unit thickness) as shown in Fig. 3. The 


(4a) 

(4b) 
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stress intensity factor is given by 

(5a) 
(5b) 


K II - 


§yb r« 


or 


?rb 


This result may be used to generate (by superposition) the stress 
intensity factor for an arbitrarily loaded crack face. 

A popular model for representing the anti-plane deformation 
field of a strike-slip plate boundary (uniform along strike) with 
a shear zone sliding under a locked seismogenic layer is shown in 
Fig. 4. The quasi -plastic shear zone is modelled as the lower 
edge crack of length a while the upper edge crack with length b 
has been used to simulate shallow creep by Tse et al (1985). Of 
course, putting b=0 is equivalent to locking the shallow crust, a 
model first employed by Turcotte & Spence (1974) in analyzing 
surface strain profiles along a line perpendicular to the plate 

boundary. The stress intensity factors at the lower and upper edge 

crack tips are given by (Tse et al, 1985). 

K - crVH 72sin(7ra/H)/(a+/3) (6a) 

and 

K m (b) - aj H y2sin(7rb/H)/(a+£) (6b) 

where a -= cos(?ra/H) and /3 — cos(7rb/H). In section 4.4 this model 
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is used in a more detailed analysis of plate boundary 
deformations . 

2.1.2 Energy Release Rate G 

The energy release rate G is defined to be the energy flux to 

the crack tip zone per unit crack length advance (per unit width 

along crack front). For mode II shearing, this definition affords 

a connection between G and (see e.g. Irwin, 1960). Referring 

to Fig. 5 and by recognizing that the process of crack extension 

is to cause the material in a small zone A £ to slip an amount 

6, = ^ — K 78(Ai-x 1 )/7r (lb) with stress reduction A a., from 

1 n 11 1 21 

K II f f 

— + a (la) to the residual friction a , the work absorption 
Jinx 

per unit crack advance and hence the energy release rate may then 
be calculated as 


G 


lim 
Ai -» 0 


1 

2Ai 



Aa 2 i( x i) 6 1 (Ai-x 1 ) dx 1 


lim 

Ai -*• o 


1_ (!-»/) 
A i fin 




J 


Ai -x. 


dx 


1 


- -13^2 K 2 

2n II 


(7) 


If other modes of deformations are involved, the additional work 
absorbed into the crack tip has to be acounted for and in that 
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For generalized plane stress deformation in mode I or II, (1-v) 
should be replaced by l/(l+v) in (8). The linear superposition is 
appropriate here because the three deformation modes are 
independent of one another directly ahead of the crack tip. This 
same form is obtained by Dmowska and Rice (1986) in an elegant 
presentation starting with an infinitesimal growth of a general 
three-dimensional crack front and considering the associated 
energy changes in the body containing such a crack. For a given 
mode of fracture, (8) explains why the Griffith's fracture 
criterion based on G is equivalent to Irwin's fracture criterion 
based on K. 


The critical energy release rate G c of the earth's crust may 
be estimated by various means, which we shall discuss in some 
detail in section 3. For now we would like to illustrate a use of 
elastic brittle crack mechanics with a simple example. Consider 
the creeping segment of the San Andreas fault in central 
California as a large crack of length 21 in an elastic plate (the 
lithosphere) under mode II generalized plane stress deformation 
(Fig 2). This is of course a crude approximation only, although 


the currently locked segments of the 1857 and 1906 ruptures 


the enhanced background seismicity near San Juan Bautista and near 
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Parkfield lend credence to such a representation. The enhanced 
seismicity and the short cycle (approximately 21 years) of 
moderate Parkfield earthquakes may be considered as local slip 
instabilities in the breakdown zones of the megascale crack. With 
an effective length St and a driving stress A a, the crack face slip 
displacement is given by (see, e.g. Muskhelishvili , 1953) 


8 


2A a 




(9) 


Equation (9) is used to generate a curve fit for field data of 
fault creep in Fig. 6, which shows four sets of fault slip rate 
measurements (after Burford and Harsh, 1980; Lisowski and 
Prescott, 1981; Schulz et al , 1982) and slip rate prediction based 
on (9). (Since we have a linear problem, the slip rate 8 is 
related to the stressing rate A a in the same manner as 8 is 
related to A a in (9)). The geodolite data set should be weighted 
more than the other sets because it reflects relative displacement 
further out (up to 5 km) on each side of the fault trace than the 
other data sets and is therefore more suitable for the two 
dimensionality of the present crack model. 

The energy release rate may be estimated from 

,2 

G - j; n(l+i /) max (10) 

Co n 


Equation (10) has been obtained by combining (4) , the plane stress 
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version of (7) and (9) with S «= <5(x-0). The field data show a 

max i 

maximum slip rate of 3.4 cm/yr near Monarch Peak. If an average 

repeat time of 160 years is assumed for great ruptures of the 1857 

or 1906 type, then ^ max ~ 5.44 m. Equation (10) then gives 
6 - 2 

G^ = 6.3 x 10 Jm for i/=0.25, /i = 35 GPa and 2 — 80 km. The 

estimated value of G may yet be higher if one considers that 

seismic ruptures occur in the shallow crust of say, 10 km in a 

50 km thick lithosphere. In that case, the (thickness -averaged) 

G c should be weighted by a factor of 5 which then results in 

G = 3.2 x 10 7 Jm' 2 . 
c 

Rice and Simons (1976) also used (10) to calculate G for a 

c 

fault creep event on the San Andreas which occurred on July 17, 

1971 and was reported by King et al (1973) . The maximum slip 

value 6 recorded by the four creep stations was 9 mm and the 
max 

length of the creep zone was reported to be 6 km. Using j/= 0.2 and 

2 -2 

M“20 GPa, Rice and Simons calculated a G - 2.6 x 10 Jm . This 

c 

is one of the lowest values for G^ reported based on field 
observations. Presumably the episodic creep events involve 
extension of slip zones in clay gouges with low strength (see 
section 3 and Tables 3 and 4) or small confining pressure at 
shallow depth, whereas seismic ruptures of the 1857 or the 1906 
types involve the breaking of both fissured and competent rocks. 

An alternative more fundamental view of the energy release 
rate, based on its basic definition but without referring to the 
linear elastic stress and deformation fields given by (1) or (2), 
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is the following: Suppose a tube of material bounded by the 
contour Y surrounding the crack tip is cut, as shown in Fig. 7. 
The energy flux going into this tube of material is used up in 
elastic straining of the material in A, in frictional work on the 
crack face L, and in supplying energy to drive the extension of 
the crack. If the crack extends at a steady speed v = da/dt, then 
for unit thickness in direction x^ , 


J T.(du./dt)dr = ^ J WdA + J o 1 


(d5 i /dt)dx 1 + G(da/dt) 


(ID 


where T - a-n is the traction vector acting on T, and W is the 


elastic strain energy density, defined as 


i: 


mm 


a . . de . . . This 


statement of energy balance presumes that any energy absorbed by 
inelastic deformation apart from frictional work, can be lumped 
into G. For steady state extension where self similarity is 
preserved (i.e., an observer riding with the moving crack tip 
always observes the same view) one can replace the time derivative 
a/at by -va/ax 1 since in that case x p = x p * vt where x^ refers 
to a coordinate system fixed in space. Also if the contour T is 
shrunk onto the crack tip, the frictional work term can be 
eliminated. Thus (11) becomes 


rio I ["!« 


T. 


au.-i 

i 


i ax 


i J 


dr 


( 12 ) 
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For dynamic crack propagation, an additional term involving 
kinetic energy dissipation must be added inside the integral 
(Cherepanov, 1968; Dmowska and Rice, 1986; Freund, 1976; Kostrov 
and Nikitin, 1970) . 

2 . 2 The J- integral 

The J- integral is defined by (Rice, 1968a) 


J 



T. 


3u. n 

l 


i 3x, 


dr 


(13) 


for two-dimensional quasi-static deformation fields in elastic 
solids. Comparison with (12) shows that G is equivalent to J in 
the limit that the contour F shrinks onto the crack tip and when 
the self - similarity condition leading to (12) holds. In general J 
may be interpreted as the excess of energy flux through T over the 
elastic strain energy absorption by the material inside T. 

The J integral has been shown to vanish on a closed contour 
containing no singularity, for an elastic material which is 
homogeneous, at least in the x^ direction (Rice , 1968a) . The 
J- integral is one component of a set of three conservation 
integrals noted by Eshelby (1957) as characterizing energetic (or 
configurational) forces on localized inhomogeneities in elastic 
solids, and its exploitation in treating crack problems was first 
pointed out by Rice (1968a, b). An energy release interpretation 
as in (12) was published by Cherepanov in 1967, (in Russian). (He 
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did not seem to realize that the J- integral was path- independent 
at that time . ) . 

For a mode II shear crack with shear traction a on the crack 
face (Fig. 8), the conservative property implies that 


J Q + j , + - Jp + •> 

v Q P P Q 


(14) 


Jq is the J- integral on a contour starting from a point Q on the 

upper face of the crack, surrounding the crack tip and ending on 

the point Q on the ■ lower side of the crack face; J is the 

Q P 

J- integral on a short contour from Q + to P + along the upper crack 


face, etc. Equation (14) may be reduced to 


Jq + J <7(36/3x^)dx^ 

Q 


f 


cr(d5/3x^)dx^ 


(15) 


where each integral already incorporates contributions from both 
the upper and lower crack faces on which n^=0 and 6 = u^ + -u^ . In 
the case where the crack faces are traction-free, i.e. 0, then 
J^Jp, resulting in the path- independent property of J. In shear 
cracks, however, a is usually non-zero. The statement (15) is 
valid for any points P and Q. If we allow the point P to approach 
the crack tip 0, the integral term on the right hand side of (15) 
vanishes and the remaining part Jp is just the energy release 
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rate G, by (12) and (13). Thus for mode II shear cracks 


G 




(16a) 


Obviously the right hand side of (16a) cannot depend on the 

f 

specific point Q. For uniform a=a , the integral term in (16a) 
is simply -o f 5q, i.e., 



(16b) 


Equation (16b) is a result of the small scale yielding condition 
for a shear crack, and affords another interpretation of J : It is 
the energy sum of crack driving force and frictional dissipation. 

Equation (16b) may be exploited to obtain estimates of the 
crack driving force by choosing contours on which the terms in 
(13) can be easily evaluated. As an illustration, Rudnicki (1980) 
estimated the G c for initiation of the 1857 California rupture 
using (16b) and the contour shown in Fig. 9. The contour is 
chosen such that the left vertical branch cuts across the San 
Andreas where it has been locked and the material there is assumed 
to deform uniformly with shear strain 7 . Similarly the point Q 
is located well inside the creep zone in central California such 
that uniform shear straining 7 ^ may again be assumed. On both of 
the vertical contours, a pure shear state assumption implies that 
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3u^/3x^=0 . (The origin and the crack tip are located 
approximately at Cholame) . The quantity 5^ represents the creep 
magnitude at point Q. The horizontal contours are chosen at 
X£ = ± h /2 where loading is essentially imposed by displacement. 
There n^=0 and 3u^/3x^=0 (due to uniform imposed displacement) so 
that no contribution is made to the J- integral. The only 
contributions come from the strain energy of the vertical 
branches. Thus 


J q - W( 7 o )h + 2 W(y r ) (-h/ 2 ) 

f 7 0 

“ h o(-y)d~t 


The fault creep 5^ may be estimated from the difference between 
the strains well inside the locked segment and well inside the 
creep zone, i.e. 


5 


Q 



7 )h 
r 


*= h 



d7 


Thus at rupture initiation the critical energy release rate, using 
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(16b) , is 


G 

c 


- h 


r 7 o r . . fi , 

] a (7) - a dy 

J v L 


2 h 



(17) 


if it is assumed that a(y) - a 1 


7-7 

r . 


for linear elastic 


behavior. It is interesting to note that (17) has the same form 

as (10), with 1 replaced by h and S by £ . The factor 7 r(l+i ^)/8 

max Q ' 

in ( 10 ) also is close to 1/2 as in (17) for v = 0.25. 

For numerical estimates, 5 ^ should reflect the part of the 

relative plate motion accommodated by the San Andreas. Minster 

and Jordan (1984) suggested a slip rate of 35 mm/yr, which 

translates to 6 ^ « 0.035 x 145m — 5.08m for a repeat time of 145 

years (Sieh, 1984) . Strain rate profiles based on geodetic 

measurements (King, N.E., personal communications, 1984) on the 

Carrizo Plain and San Luis net decay from the fault trace and 

appear to flatten out at approximately 60 km. Thus 60 km seems to 

be an appropriate value for h. Using a crustal averaged shear 

6-9 

modulus fi «* 35 GPa, (17) gives G c to be 7.5 x 10 Jm" . Again 

7 - 2 

this value would be increased to 3.8 x 10 Jm if we consider 
that seismic ruptures occur to a depth one-fifth the lithospheric 
thickness. These values are slightly larger than Rudnicki's 
original estimate because of different values assumed for S and 
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6 7 - 2 

for /i . They are consistent with the 6 . 3 x 10 - 3 . 2 x 10 Jm 

estimate using surface creep rate data from central California 
(section 2.1.2). The calculated value of G should be considered 
as an order of magnitude estimate only because of simplifying 
assumptions. For example the model shown in Fig. 9 assumes that 
the creep zone extends indefinitely north of Cholame. This can 
only be an approximation of the finite length 160 km) of the 
creeping section of the San Andreas. 

This section reviews shear fracture mechanics. The fracture 
parameters K, G and J are introduced, and their relation to one 
another emphasized. Applications of fracture mechanics to model 
plate boundary deformation and various techniques in extracting 
fracture parameters based on seismic and geodetic observations are 
demonstrated. Elements of this section form the basis for further 
discussions of the mechanics of shear rupture in sections III and 
IV. 

III. SLIP -WEAKENING MODEL OF SHEAR RUPTURE 

Laboratory observations from tri-axial tests in rocks 
indicate a complex breakdown process in the localized shear band 
in the post-peak regime. This breakdown process may involve 
buckling of slender columns in grains segmented by microcrack 
arrays, kinking in plate -shaped grains and rotation and crushingof 
joint blocks as seen under SEM (Evans and Wang, 1985). On a 
larger scale, direct shear testing of rock joints indicates 
shearing off and crushing of asperities in jointed rock mass, the 
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micromechanics being sensitive to the normal stress applied across 
the joint (see, e.g., Coulson, 1972). Both tri-axial rock 
specimens and direct shear jointed rock specimens in the 
laboratory show a decreasing shear load carrying capacity as a 
function of the amount of sliding. Some samples of such 
experimentally measured slip-weakening curves for initially intact 
rock specimens are shown in Fig. 10. Slip-weakening curves for 
jointed rock specimens are shown in Fig. 11. These slip -weakening 
o-8 relations define simple constitutive laws governing shear slip 
behavior. While direct use of experimental results in the field 
for earthquake faults may not be appropriate because of 
differences in size scales of fault zone structures in comparison 
to laboratory scale specimens, one may expect that certain 
behavior of fault slip may be governed by similar slip -weakening 
relations. This is so because for a fault to exhibit seismicity, 
its strength must degrade with on-going slip. 

In the following discussion (Section 3.1) a general 
constitutive model for the slip-weakening process is introduced. 
The model is an extension to shear faulting by Palmer and Rice 
(1973) and Ida (1972) of the well-known cohesive zone models of 
tensile fractures developed by Barenblatt (1962) and Dugdale 
(1960) for metal. (The cohesive zone model was also used to 
describe crack extensions in concrete by Hillerborg et al (1976) 
and Li and Liang (1986). In this case the breakdown process 
involves the joining of discontinuous microcracks and pull-out of 
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the aggregates from the cement matrix) . The slip -weakening model 
will then be applied to an analysis of stability of a simple 
sliding system (Section 3.2). Relationship between the 

slip -weakening model, the J- integral and the elastic brittle crack 
model will then be described in Section 3.3. In Section 3.4 we 
shall discuss some estimates of the parametric values in the 
slip-weakening a-S relations from laboratory tests of rocks, rock 
joints and overconsolidated clay and from in-situ field 
observations of natural fault behaviors. 

3 . 1 Slip-Weakening Constitutive Model 

A simple general form of the si ip -weakening constitutive 
model may be written as 


a - f (5 , a; , T) (18) 

where o' ** a - p is the effective normal compressive stress 
n n 

(normal stress reduced by pore pressure p) acting across the 
slip surface, and T is temperature. A schematic plot of (18) is 
shown in Fig. 12 which indicates a continuous decay of strength 

from a** to a at large slip beyond a critical slip displacement 

6 *. 

P f 

Triaxial tests in rocks by Wong (1986) suggest that a and a 

P f 

both increase in a manner such that a - a first increases but 
then decreases with increasing , indicating a transition from 
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brittle to ductile deformation. In a separate series of tests at 

constant normal stress, Wong (1982b) shows that the stress drop 

P f P 

a -a decreases with temperature. The ^ and T dependence of a 

f 

and a is illustrated in Fig. 12b, c. These considerations of 
normal stress and temperature dependence are important when the 
constitutive law is applied to the earth's crust, as in the work 
of Stuart & Mavko (1979), Li and Rice (1983a, b) and Li and Fares 
(1986) . It is assumed in the slip -weakening relation that 
unloading and reloading from the weakening branches occur along 
vertical paths (Fig. 12a), i.e. no reverse sliding accompanies a 
load removal . 

The model as described has no dependence on slip rate, 
although recent rock experiments by Dieterich (1978,1979), Ruina 
(1983) and Tullis and Weeks (1986) indicate that frictional 
sliding behavior has slip rate and state dependence. An 
implication of ignoring rate and state dependence is that slip 
events could not be repeated on the same surface since no 
mechanism for restrengthening is available. Thus the 

si ip -weakening model is unable to describe transitions from one 
earthquake cycle to another. Even so, the model is capable of 
simulating sliding behavior quite adequately in most situations. 
It is in fact a more general description of shear slip phenomenon 
than the elastic brittle crack model, which we shall reveal as a 
limiting case of the slip -weakening model in Section 3.3.1. 
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3 . 2 Stability Analysis of a Single-degree-of-freedom System 
3.2.1 Sliding of a spring-block system 

An earthquake may be regarded as a loss of stability in slip 
on the fault surface. If the slip weakening a-S relation is a 
fair constitutive description of the fault zone, then it should, 
in the minimum, allow a sliding surface to initiate a rapid slip 
event subsequent to some stable slip. It is important to note, 
however, that whether an unstable event is generated or not 

depends not only on the constitutive relation of the sliding 
surface but also on the stiffness of the system (loading system 
and the body which contains the sliding surface) through which 
load is transmitted to the surface. To see this more explicitly, 
consider the mechanical behavior of a simple spring-block system, 

as shown in Fig. 13a. The block is assumed to be rigid and the 

sliding surface is governed by a slip-weakening relation (solid 
line) shown in Fig. 13b, i.e., the shear stress a acting on the 
sliding surface and the block movement 6 follow such a 

relationship. The block is loaded through a spring which is 
pulled forward by the amount 6 . The normal stress acting on the 
block, as well as the temperature on the sliding surface, are 
assumed to remain constant during the sliding process. 

The force equilibrium equation governing the system can be 
written as 

T - a (19) 

where T is the spring force. The load and load point displacement 
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are related by 


.T ■= k (6 -6) . (20) 

This incorporates information of the structural stiffness of the 
medium (the spring) through which loading is applied to the 
sliding surface. Equations (19) and (20) may be combined and 
rewritten as 


a = -k£ + k«S (21) 
o 

which then defines an unloading line of the loading system with a 
negative slope -k in the a- 8 space (Fig. 13b), and with an 
intercept kS Q on the vertical a-axis. Equation (21) expresses 
that equilibrium of the system is satisfied on any point of this 
unloading line. However, for each unloading line shown, only its 
intercept with the o-8 curve can be the true equilibrium point 
since the sliding surface is governed by the constitutive relation 
a « o(8) , i.e. , 


o(8) - -k5 + k5 Q (22) 

Thus a series of equilibrium points A, B, C, D may be traced as 6 
is increased by pulling on the spring. The unloading lines 


illustrated have been drawn at equal (vertical) intervals such 
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that S q is increased uniformly in time simulating a steady load 
point displacement. The block displacement (5 -»£ -+$ ...), 

however, accelerates with time. For the o-S relation and spring 
stiffness k shown in Fig. 13b, equilibrium could be maintained 
only up to Point E. In the sequence from A through E the force in 
the system rises to a peak (at D) and then decreases. At point E 
the system becomes unstable in the sense that an infinitesimal 
increase in S causes a sudden jump in S (the block shoots 
forward from 6 ) accompanied by a stress drop. This is 
schematically illustrated in Fig. 13c. The final equilibrium 
position can be any one of the states F., G or H. These points are 
constrained by the fact that unloading of the spring must follow 
the unloading line EE at instability, and that the sliding 
surface must still obey the slip -weakening law (including the 
rigid unloading branches, Fig. 12a). Furthermore, the energy loss 
from the spring must be converted into work of sliding the 
surface, which implies that 


1 

2 


( “E + V 




a(<$)dS 


(23) 


which is how the point H is defined. However, if energy - is 
partially lost through seismic radiation (or heat generation other 
than that related to the frictional work on the right side in 
(23)), then the final resting equilibrium position may be at F or 
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G, and the equal sign in (23) should be replaced by > when 5 is 
replaced by 5p or 5 As mentioned earlier, the slip weakening 
constitutive framework do not allow a natural restrengthening of 
the slipped surface. Hence reloading from G or H would follow 
first the rigid branches and then along the residual friction 
plateau. In such a system, a single instability is allowed but no 
repeated events can be generated. 

The question of how much stable sliding and at what load 
level instability sets in could be addressed by considering 
springs of different stiffnesses. The case of an infinitely stiff 
spring (vertical unloading lines) is analogous to applying the 
load directly onto the rigid block such that the load point and 
the block move stably together (6=6^). In the case of an 
extremely compliant spring (k-*o) , the unloading lines are close to 
horizontal and the instability point would be at the peak of the 
si ip -weakening curve. In the case of a stiff spring whose 
stiffness is greater than the negative of the slope at any point 
of the weakening branch, as shown in Fig. 13d, no instability 
would occur. The stress and slip development are schematically 
shown in Fig. 13e for this case. Note that while the slip 
accelerates, its rate does not approach infinity, as in the case 
of a true instability (Fig. 13c). This may be a useful conceptual 
model for what is known as "slow earthquakes" (see, e.g. Sacks et 


al, 1978). 

The considerations of where on the weakening branch 
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stability is lost explain why it is necessary to use very stiff 
machines in the laboratory (the loading assemblies act as part of 
the springs ) if one is interested in tracing the weakening 
branch in a direct shear test. Of course the remarks here apply 
to tension or compression loaded -specimens just as well (see, e.g. 
Jaeger and Cook, 1969). 

Often, in real geophysical systems, the load transmitting 
medium behaves viscoelastically . A more detailed discussion of 
such a medium will be given in section 4.2. For now we consider a 
simple single-degree-of-freedom system, and incorporate the 
viscoelastic behavior in the form of a standard linear element as 
shown in Fig. 14a. This element has thfe property that for long 
time response, or under d (5 S) /d t » 0 relaxed condition, the 
system stiffness is given by k^ . For short time response, or 
under 3 ( J 5^-5 | )/3t-+» unrelaxed condition, the increased stiffness 
approaches k^ + . In between these two limits, the dashpot 
modulates the contribution of k^ to 'the total stiffness. In 
reference to Fig. 14b the system undergoes stable deformation up 
to the point I, where the unloading line associated with a relaxed 
stiffness k^ is tangent to the slip -weakening curves. Once rapid 
acceleration is initiated, the dashpot is activated and 
effectively stiffens the loading medium. At this point the system 
is self driven in the sense that the block continues to slide to 
the right even if the load point has stopped moving, but the 
motion may still be considered quasi-static, following the path I 
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to D. At point D where the element stiffness has reached its 
maximum limit k^ + dynamic instability sets in, an analogue of 
an earthquake . The period of deformation corresponding to I to D 
may represent a precursory period when anomalous activities 
associated with rapid straining are revealed. The viscoelastic 
element sets the time scale of this precursory period. 

The discussion presented in this sub-section provides a 
general conceptual framework for stability analysis of any 
slip -weakening or strain-softening system. The instability 
delayed mechanism could be associated with viscoelastic 
stiffening as described above, or it could be associated with 
drainage responses in a poro- elastic medium. 

3.2.2 Application to Fault Stability Analysis 

As an illustration of some of the saliant points raised in 
the above instability analysis, we digress from the spring- dashpot 
si ip -weakening model to consider a more realistic time -dependent 
fault system. Li and Rice (1983a, b) analyzed the stability of 
stressing of a seismic gap zone in which progressive failure 
eventually lead to an earthquake at a strike-slip tectonic plate 
margin. The actual problem involves a seismic gap zone of length 
22 in an elastic lithospheric plate underlain by a viscoelastic 
foundation. The lithosphere is assumed to undergd plane -stress 
deformation and is coupled to the asthenosphere in the form of a 
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simple foundation (Fig. 15a). A modified Elsasser model (Rice, 
1980, Lehner et al, 1981) is used to include the resistance (to 
rapid deformation) due to viscoelastic coupling as body forces in 
the lithosphere. Further, the driving stress a°(t)-a(t) (averaged 
over the lithospheric thickness, Fig. 15b) is assumed to be 
uniform over the gap zone. This causes the thickness -averaged 
slip displacement to be distributed as in (9) for a crack model. 
Li and Rice (1983a, b) related the average of this slip over the 
seismic gap to the driving stress using the representation 

S(t > - f 

which is an extension of (21) to incorporate the viscoelastic 
effects of the stress transmitting medium. Here the driving 
displacement kS Q has been replaced by the driving load term a° and 
of course k is the inverse of the compliance C(t) in (24) . Indeed 
this correspondence may be easily seen in the long-time (relaxed) 
elastic limit, in which case (24) becomes 

S(t) - C(«) ja°(t) - a( t)j (25) 
and in the short time limit, in which case (24) becomes 


C(t-t ) 


dt 


(t ) - cr(t ) dt 


(24) 


d5 - -C(0)do 


(26) 
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To further relate a to the physical driving mechanism - the 
relative plate velocity - we may consider the time derivative 
of (25). In that case, the slip rate d6/dt averaged over several 
earthquake cycles must correspond to and the stressing rate 
da / dt must average out to zero. Thus do°/dt -= V /C(«>), consistent 
with the loading term kS Q in (21) . This direct interpretation of 
a° in terms of V ^ was noted by Tse and Rice (1986) . 

The compliance function C(t) is shown in Fig. 16 for 22 ** H 
and 22 « 5H, obtained from numerical inversion of C from the 
Laplace transform space (Li and Rice, 1983a, Appendix B) . The 
zero time and infinite time limits of C(t), however, are derivable 
in analytic forms 


C(0) 


H 

(l+i v)£p 


I. 


22/Ur 


erf 


y?]‘ 


d 4> 


(27a) 


C(«) 


8L 

jt(1+i 


(27b) 


Naturally, the compliance at infinite time corresponds to a plate 

with a completely relaxed foundation and could be obtained 

directly from (9) by taking the average over -2 to +i in x^. 

2 

(However, (27b) is actually 16/n times the exact result from (9), 
due to an approximation in representing a finite length crack by a 
semi- inf inite crack with a uniformly loaded finite portion from 
the crack tip (Lehner et al , 1981). A way to compensate for this 
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discrepancy is to regard i as an effective crack length equal to 
2 

n /16 of the actual length) . The compliance at time zero 
corresponds to a fully coupled lithosphere/asthenosphere system so 
that C(®)>C(0) . In addition, the compliance may be expected to 
increase monotonically with seismic gap zone length i and to 
decrease with shear stiffness fi. 

As explained earlier in connection with the standard element, 
the effect of an increasing stiffness (or decreasing compliance) 
with increasing slip velocity is to delay the final instability. 
Li and Rice (1983a, b) analyzed the details of this precursory 
stage by solving (24) together with a crustal scale (averaged over 
the lithospheric thickness) a-8 relation of plate boundary 
deformation. It should be noted that this a-8 relation is not a 
material constitutive law as for the si ip -weakening model, even 
though it exhibits similar behavior of decreasing a with 
increasing 8 as the slip zone penetrates into the seismogenic 
zone, as described below. 

The crustal scale a-8 relation is derived based on an 
anti-plane strain analysis of an edge crack strip (Fig. 13c or 
Fig. 4 with b=0) representing the deformation behavior averaged 
over the seismic gap zone. In this case a and 6 are related 
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parametrically through the crack length a 


/iG c (a)/H tan(jra/2H) 


1/2 


(28a) 


S = (4HA/i) ^G c (a)/H tan(?ra/2H)j ^ In Fl/cos (wa/2H) j (28b) 


by requiring that a be of magnitude that just meets the fracture 

criteria at crack depth a. Thus G^Ca) is a prescribed quantity 

which should reflect the changing fracture property of the shear 

zone in the Earth. Section 3.4.1 describes experimental 

observations of the effect on G ^ due to changes in temperature and 

pressure. Li and Rice (1983a, b) choose a Gaussian distribution 

with depth (Fig 17a) , with parameters adjusted to fit typical 

focal depths and the thickness of the seismogenic layer. Thus, 

for example, the maximum of G , or G , lies at the seismogenic 

c max 

depth . 

It may be noted that (28a) is consistent with the stress 

intensity factor calculation (6a) and (8) for an elastic brittle 

crack model. Indeed (28a) affords an estimation of the critical 

energy release rate in the earth's seismogenic zone. Based on an 

average stress drop of 30 bars reported to be typical of great 

plate boundary ruptures (Kanamori and Anderson (1975)) and a 

focal depth of about 10 km appropriate for great ruptures on the 

6 2 

San Andreas, G was constrained to 4 x 10 J/m . Li and Rice 
max 
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(1983a) also suggested that this choice of G c produced slip 
magnitudes of 2.5 - 4.5m consistent with observed seismic slip 
magnitudes of great California ruptures. 

The resulting a-S relation based on (28) and on the Gaussian 

distribution of G c is shown in Fig. 17b. Stability analysis of 

the plate boundary may follow the graphical analysis of the 

single -degree -of- freedom system shown in Fig. 13 or 14, at least 

up to the point of initial instability. This corresponds to the 

state when the compliance C(t) just drops below that of C(»), 

after which the numerical solution of (24) becomes necessary. 

Naturally (24) must now be regarded as an integral equation for 6, 

when <7 inside the integral is expressed in terms of 6 through 

(28) . Li and Rice gave numerical solutions following through the 

process from peak stress, through initial and dynamic instability. 

The time -evolution of a, S, and a are shown in Fig. 18, for two 

different loading rates defined by the parameter R ** t <7°/(K m /,/H) , 

in which t is the relaxation time of the asthenosphere (see e.g. 

Lehner et al , 1981 or Li and Rice, 1983a) and K is the fracture 

m 

toughness corresponding to G » It may be noted that while the 
plate boundary stress a is decreasing, (Fig. 18b) the average a 
and 5 (Fig. 18a, c) are increasing and their time derivative 
reaches infinity at dynamic instability. The time scale of the 
self -driven progression from state I to D depends, among other 
things, on the relaxation time t of the asthenosphere. 
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While the model described above gives a plausible 
representation of the gross phenomenon of the precursory processes 
leading to loss of stability, i.e., an earthquake, the 
single-degree-of-freedom model is obviously an over-simplification 
of the source mechanism. The earthquake source process may 
involve a local nucleation followed by subsequent spreading of the 
fault surface along strike. This is even more likely when one 
considers the possibility of spatial variation of material 
properties and geometric features so that a, 5 and a may be highly 
location dependent on the fault. To include along- strike 
variation, Tse et al (1985) and Li and Fares (1986) analyzed a 
multiple -degree -of -freedom system to be described in section 4.4. 
Another inadequacy is related to the assumption of the elastic 
brittle edge crack model to represent the slip penetration which 
leads to the unavoidable necessity of the loss of stability as the 
slip zone approaches the ground surface; i.e. no stable continuous 
creep could be simulated. We shall reexamine this issue after 
introducing the relationship between the crack model and the slip 
weakening model in the following section. 

3 . 3 Slip-Weakening Model. J-integral and Elastic Brittle Crack 
Model 

The behavior of shear failure as represented by a 
single-degree-of-freedom system detailed in the section above is 
perhaps suitable for small surfaces such as in typical laboratory 
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specimens. However, such representation is plainly inaccurate for 

a large surface such as a natural fault. This is because the 

slippage at each position on a fault surface may be quite 

different: At some position which has undergone extensive 

f 

sliding, the stress level may be close to a , whereas other 

p 

positions may still be high up close to a on the slip -weakening 

curve. This distributed slip situation contributes to the 

phenomenon of stress concentration, especially at material points 

where slip has barely begun (e.g. on the positive sloping branch 

of the slip -weakening curve, Fig. 12a). Indeed the business of 

fracture mechanics has, to a large extent, to deal with the 

intensity of such stress concentrations. Thus the elastic brittle 

crack model may be considered an extreme member of the 

slip-weakening model: Outside the crack the material remains 

elastic, whereas inside the crack the material has all slid down 

f 

to the residual friction level a , with the exception of a small 
zone at the crack tip. The smallness of this zone corresponds to 
the so called small scale yielding condition in elastic brittle 
crack mechanics. This condition may be addressed within the 
context of the slip-weakening model if we assume that the 
inelastic deformation within this small zone is governed by the 
slip-weakening constitutive relation. 
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3.3.1 Relationship Between Energy Release Rate G and Slip- 
Weakening Model Parameters 

We show in Fig. 19a the stress distribution near the tip of a 

mode II sliding surface. The zone of size w contains the stress 

P f 

degradation from peak strength a to residual friction a , at 
which slip 8 has reached the critical value 8 . The weakening 
branch of the si ip -weakening curve is shown in Fig. 19b. 

We now apply the J-integral with a contour going along the 
lower and upper crack faces surrounding the breakdown zone, 
beginning and ending at a point Q located well beyond x^ ■= -w 
(Fig. 19a). Recognizing that no crack tip singularity exists due 
to the presence of the breakdown zone, Eq. (16a) implies 


0 


J Q + 



o(35/3x^)dx^ 


(29) 


where a = a(S) according to the slip -weakening relation. The 
integral term may be rewritten as' 



(30) 


by realizing that a is a single valued function of S (if no 

unloading occurs) and hence (35/dx^)dx^ *= d5, and that the square 

* 

bracket on the right hand side of (30) vanishes for 8 > S . Thus 
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(29) and (30) leads to 



( 31 ) 


which may be interpreted as follows: The excess energy flux 
(above doing frictional work) made available balances the energy 
absorption in the breakdown process for slip zone extension. 
Equation (31), with the left hand side interpreted as a crack 
driving force and the right hand side interpreted as a fracture 


resistance , then affords a criterion for propagation of the slip 
zone. It is useful to note that the quantity (J^ - cr^S^) cannot 
depend on the particular point Q (since the right hand side of 
(31) is independent of the point Q) , which suggests that 

f 

(Jq - a 6^) is a path- independent parameter for cracks with crack 
face tractions, with the stipulation that the point Q be outside 


the breakdown zone. Palmer and Rice (1973) used (31) to evaluate 


the criterion for the propagation of a shear band in 
over-consolidated clay in a long shear box and for the extension 
of a slip surface in a soil slope loaded by gravity in response to 
a step cut in the slope. Indeed Rudnicki's calculation of for 
initiation of the 1857 rupture on the San Andreas described in 
section 2 is an extension of the shear box analysis by Palmer and 
Rice . 


It should be noted that the derivation of (31) makes no 
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assumptions about the size of the breakdown zone. To make contact 
with elastic brittle crack models, however, (16b) based on an 
infinitesimal breakdown zone (consistent with the uniform a — 
assumption) may be combined with (31) to give 

G - J q o(S) - o f j d6 (32a) 

within the context of the slip-weakening model. This integral is 
just the shaded area under the slip-weakening curve (Fig. 19b) and 
may be rewritten as 



(32b) 


in which the nominal slip distance 6 is defined as 


S 


1 


P f 
a -a 



f 

a 


d 6 


Thus in the limit when the size of w is small, the slip -weakening 

model is consistent with the elastic brittle crack model with 
P f — 

G c ** (o - a )6. At incipient faulting, the critical energy 
release rate G is now interpretable in terms of the product of 
the stress drop and the nominal slip, parameters which describe 
the crack tip breakdown process. 
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3.3.2 Estimation of the Breakdown Zone Size w: 

The breakdown zone size w at crack initiation may be 

estimated from the fact that the net stress intensity K must 

J net 

vanish for no stress singularity at 0, i.e., 

K — K . . , + K, , - 0 (33) 

net applied Ddz 

where K .^ e{ j i s th e stress intensity factor due to applied load 

which is just equal to the critical value K c at initiation of 

crack growth, and - is the stress intensity factor due to 

excess (over friction) shear resistance in the breakdown zone. 
Clearly, is a negative quantity and may be calculated if the 

stress distribution in the breakdown zone is known. In general 
this information requires the solution of a singular integral 
equation which we shall describe in section 4.1. As an estimate, 
Palmer and Rice (1973) used an inverse method in which a linear 
stress distribution is assumed, and the resulting a-S relation is 
shown to resemble an actual si ip -weakening curve (Fig. 20). In 
this case 


o(r) 


(/ - a f ) (1 - r/w) 


(34) 



48 


and using (5b) and superposition, 



g(r)-g f 


dr 


Thus (33) together with (34) and (35) implies 



[ 


K 

c 


P f 
a -o 




2 


or using (7) and (32b) 


(35) 


(36) 


w 


97T(l+t/) 

16 


pt6 


t P fv 
(a -a ) 


(37) 


Full numerical solution of the singular integral equation 

mentioned earlier indicates that (36) gives a good estimate within 

10% error if the slip-weakening curve has a linear decaying shape 

(Li and Liang, 1986). However, (36) is inadequate for a material 

with an exponentially decaying slip-weakening curve with a long 

★ 

tail (i.e. large 6 ). 

This section summarized the connection of the slip-weakening 
model with the J- integral in general, and with the elastic brittle 
crack model in the limit when the breakdown zone size w is small. 
Indeed in the context of the si ip -weakening model, w must be 
smaller than all other characteristic dimensions (crack length, 



49 


distance of crack tip to boundary, etc.) in the geometry of the 
body in order to satisfy the small scale yielding condition 
required in the use of the simple elastic brittle crack model. 
Furthermore, since this length is derived from material properties 
(37), it creates a problem in geometry scaling in laboratory model 
studies (e.g. centrifuge ’studies), a difficulty first noted by 
Palmer and Rice (1973) . 

It should be clear now that for most laboratory size 
specimens, w is generally relatively large and may well exceed the 
dimensions of the slip surface. In this case the 
single-degree-of-freedom system (Section 3.2) gives a good 
description of the slip behavior. In the field, linear elastic 
fracture mechanics may be used whenever w is small enough. In the 
following section, we examine the applicability of elastic brittle 
crack mechanics to a plate -boundary model in light of the 
discussions presented above. 

3.3.3 Crustal Scale Applications 

We now take up the question of whether the aseismic shear 
zone below the seismogenic layer at a plate boundary may be 
modelled by an anti-plane strain mode III elastic brittle edge 
crack, as shown in Fig. 4 (with b«0) , assuming that we accept the 
shear zone as indeed confined to a narrow width even as it 
approaches the base of the lithosphere. This question may be 
addressed in the context of a slip-weakening model of the shear 
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zone proposed by Stuart (1979a ,b) relating the local excess (over 
friction) shear strength a to the local anti-plane slip 
displacement 5(^U2‘ U 3) i- n t ^ ie form 

a *= o(z , S') — S exp 

Equation (38) has the same Gaussian variation with depth z as 

assumed in the crack model (Fig. 17a) with a reaching peak value 

at z q and a spread measured by d. Furthermore in connection with 

the slip -weakening terminology introduced earlier, S is the 
P f 

strength drop ( a -a ) and A is a measure of the critical slip 
displacement 6 . To make further contact with the crack 
model, (32a) requires 


2 2 
(z-z o ) Z /d Z 


exp 




2 2 
8 / A 


(38) 


G (a) - [ a(z-H-a, S)dS (39) 

C **o 


and the maximum value of occurring in the seismogenic zone is 
given by 


G 

max 


£ 


exp[-5 2 /A 2 ] d5 - ./ir SA/2 


(40) 


The size of the breakdown zone w may be estimated by (36) adapted 
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to the anti -plane mode, i.e. 


9tt ^max 

16 s 2 


( 41 ) 


For an estimate, supposing a strength drop S of 500 bars, and 

6 - 2 

using the previously estimated G ■= 4x10 Jm , the critical slip 
distance A is calculated from (40) to be A = 90 mm and the 
breakdown zone size w « 100m for n — 35 GPa. This breakdown zone 
size is much smaller than a, H-a and d which are generally greater 
than several kilometers such that the use of elastic brittle crack 
model would seem to be justified. While the strength drop must be 
greater than earthquake stress drop values averaged over the whole 
fault plane and its value of 500 bars appears to be consistent 
with the various other estimates (see Table 4), a lower value of 
say S - 100 bars (Aki's (1979) estimate for the 1966 Parkfield 
earthquake) would make A = 450 mm and w » 2 . 5 km which is still 
smaller than the characteristic dimensions in the problem but 
definitely approaching the limit of validity of the elastic 
brittle crack model. In this case it may be more suitable to 
carry out the analysis using the slip-weakening model (38). (See, 
e.g., Stuart (1979a, b) and Stuart and Mavko (1979)). In the last 
reference, the authors found that stable sliding can be attained 
by increasing the critical slip distance A and hence the size of 
the breakdown zone (see, e.g. Fig. 4a in Stuart and Mavko (1979), 
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with large A corresponding to the lower right hand corner) . The 
stable sliding phenomenon cannot be simulated by an elastic 
brittle crack model for the plate boundary. 


3 .4 Laboratory and Field Estimates of Shear Fracture Parameters 

3.4.1 Laboratory estimates of shear fracture parameters 

Rice (1980) showed that laboratory triaxial test data on 
rocks could be used to deduce shear fracture parameters. Consider 
the specimen with a throughgoing fault (pre-existing saw-cut or 
post-peak localization of shear deformation) at an angle (V2-0) 
to the major loading axis a ^ (Fig. 21a). The specimen deforms 
under confining pressure . During the test, a o ^ and the 
axial shortening AL is continuously monitored. The resulting 
curve (a^-<7 ^) vs. axial shortening AL including the softening 
branch (Fig. 21b) must be stably measured. A stiff machine (or a 
cyclic technique as used by Wong (1982a, b)) is required to prevent 
instability. The stability analysis of a single-degree-of-freedom 
system described in section 3.2 is applicable to such laboratory 
tests since the breakdown zone size w is generally much larger 
than the dimensions of the sliding surface. For example, Rice 
(1980, 1984) computed w for a series of triaxial tests conducted 
by Rummel et al (1978) in the range of 0.8 -1.2m (corrected to 
constant normal stress) for initially intact specimens. Most 
laboratory specimens have dimensions much smaller than this size. 
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Thus it may safely be assumed that sliding occurs simultaneously 
at every point of the fault of the specimen. The method of 
obtaining the slip -weakening curve from the laboratory data is 
shown graphically in Fig. 21b, d. The shear stress c and slip 6 
may be computed from 


a r a 3 


sin 26 


(42a) 


6 = AL /sin 
s 


(42b) 


where AL^ is the relative displacement of the sliding surface in 

the axial direction (Fig. 21c). However, the values of peak stress 
P f 

a and the residual stress a are affected by the normal stress 
acting across the fault, and the normal stress changes during the 
test according to 


a 2 ■ "l + "3 

a ■= r cos 26 + ^ 

n 2 2 


Rice (1984) suggested a correction procedure based on the Mohr 
circle diagram, A similar approximate scheme for reducing the raw 
data to that corresponding to constant normal stress was detailed 
by Wong (1986) , who found that the constant stress correction 


reduces the uncorrected value of G by approximately a factor of 
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Based on the above technique, Wong (1982a, 1982b, 1986) 

calculated fracture parameters from several series of tests. His 

test results (and some from other investigators) are summarized in 

Table la. Wong's studies indicate a decreasing trend in the 
P f 

strength drop a -o with increasing temperature, suggesting a 

transition from brittle to ductile deformation. Fig. 22 shows a 

composite of two series of tests, one for San Marcos gabbro and 

the other one for Fichtelbirge granite conducted at constant 

temperature. G appears to first increase with confining stress 

P f 

up to 0.55 GPa and then decrease. The G and (a -a ) variations 

c 

with temperature and normal stress are consistent with the 
observation of seismicity confinement in the shallow crustal layer 
below which quasi-plastic behavior dominates. Conducted at 
crustal scale confining pressures, these data, while still 
incomplete, are perhaps the first experimental qualitative 
evidence in support of the depth variation of slip-weakening 
parameters assigned in fault stability analysis by Stuart 
(1979a, b), Li and Rice (1983a, b) and Li and Fares (1986). See 
section 3.2.2, equation (18) and section 3.3.3, equation (38) for 
more details. The data summarized by Tse and Rice (1986), while 
phrased in terms of instantaneous and long term rate 

sensitivities, suggest a similar variation of strength drop 
potential . 

Apart from triaxial test results, slip-weakening fracture 
parameters have also been reported by Okubo and Dieterich (1984) 
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based on bi-axial tests of large scale 2m long simulated faults. 

At a normal stress of 0.6-4 MPa, they found that G c ranges from 

0.1-2. 4 Jm for Sierra white granite with prepared surface 

* 

roughness of 0.2/xm and 80/im. The critical slip displacement 5 is 

reported to be in the 0.2-10/im range for the smooth samples and 

— * 

8-40jjm in the rough samples. The values for and 5 (« 0.55 ), 

P f 

as well as the strength drop a -a , (Table lb) , are several orders 

of magnitude smaller than the corresponding values from triaxial 

P f 

tests. This lower G c and (cr -a ) may partly be due to the lower 
confining pressure under which the biaxial experiments have been 

k 

carried out, but the smaller 6 is probably related to the reduced 

surface roughness. Okubo & Dieterich (1984) reported no 

dependence of 5 on normal stress. 

Direct shear tests on rock joints have been conducted by a 

number of investigators (e.g. Coulson, 1972; Goodman, 1970, 

Barton, 1972) . Many of these studies have focused mainly on the 

P f 

effect of normal stress on a and a , with little information on 
the critical energy release rate or the critical slip displacement 
reported in the literature. However, Yip (1979) collected data 
from many rock joint tests and found a wide variation in the value 
of 5, with an averge of 0.9 mm. Like Okubo and Dieterich (1984), 

— k 

he found that this average value of 5 (or 5 ) does not appear to 
depend on normal stress. While this value of 5 is substantially 
larger than that for intact or sawcut rocks, the critical energy 
release rates G c from these rock joint experiments (based on 
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8 - 0.9mm) are lower than that for intact rocks. This is 

P f 

presumably partially due to the low stress drop (a a ) 

associated with low applied normal stress in the direct shear 
tests. Slip -weakening model parameters for some rock joint tests 
are summarized in Table 2. 

Table 3 contains si ip -weakening model parameters for 

over-consolidated clay from direct shear tests. While the stress 

drops (ct^ - a are relatively small, the 8 values are quite 

large, on the order of several mm, as opposed to less than 1 mm 

for both rocks and rock joints. This observation may have 

important implications for natural fault behavior when clay gouges 

are involved in the slip -weakening process. 

3.4.2 Field Estimates of Shear Fracture Parameters 

We have already discussed three methods of estimating in-situ 

slip -weakening model and fracture parameters from geodetic 

observations. The first method is based on the elastic brittle 

crack model and creep rate data from the San Andreas in central 

California (Section 2.1.2). There we obtain estimates for G of 

c 

6-2 7 - 2 

6.3 x 10 Jm to 3.2 x 10 Jm . The second method based on the 

6 - 2 

J- integral analysis gave an estimate of = 7.5 x 10 Jm to 
7 -2 

3.8 x 10 Jm for the 1857 Ft. Tejon rupture in California 

(Section 2.2). The third method is based on the anti-plane strain 

edge cracked strip model of Li and Rice (Section 3.2.2) which 

6 - 2 

gives an estimate of G £ - 4 x 10 Jm . Here we shall introduce 
another rather general technique for estimating fracture 
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parameters in the slip -weakening model. In this technique, 
seismological and other earthquake parameters such as rise time 

A 

t , rupture length 22 and average fault slip S are needed. Of 
course the average slip may be obtained from the seismic moment 
and fault surface area (see, e.g. Aki and Richards, 1980). Still 
other techniques and estimates for natural faults are summarized 
in Table 4. 

Consider a strike- slip earthquake rupture of length 22 (at 
least several times the plate thickness) as a mode II shear crack 
(Fig. 2). The average slip along the length of the rupture 

--2<x^<2 may be obtained by integrating (the plane strain version 
of) (9) to get 

a 

* - h 5(x i )dx i - J f 1 G (44) 

in which the relations of the stress intensity factor to stress 

drop (4b) and to energy release rate (7) have been used. The 

P f — 

product {a -a )S in the slip-weakening model may then be 

A 

calculated using (32b) and (44) once 6 is determined from 
seismological observations. As an additional constraint, 

information on the rise time t of an earthquake rupture may be 

used. The slip at each point of the fault (- 2,2 ) at time t may be 
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assumed to have the form 


S( t) = 


(a F -a f )i 


1 



(45) 


during rupture, consistent with Brune's (1970) fault model. The 
coefficient in front of the square bracket in (45) may be obtained 
from dimensional considerations and is the shear velocity. 
Assume further that the slip displacement reaches the average 

A 

fault slip S at t-t , i.e. 


S(t r ) - 6 


(46) 


combining (45) and (46) gives the strength drop in terms of the 
average slip 


P f ftS 

° - ° " ~1 

and the critical slip 6 may be computed from (32b) after solving 

P f - P f 

for G from (44) and a -a from (47). Once S and a -a are known, 

the breakdown zone size w may be calculated using (37) . As an 

illustration, for the 1976 Turkish earthquake, Purcaru and 

Berckhemer (1982) gave the following earthquake parameters: 

A 

2 H - 55 km, 6 - 2.45m, t ■* 1.5s. For n - 35 GPa, v - 0.25 and 

6 - 2 

y9 = 3 . 5 km/s , these earthquake data lead to G^ — 6.5 x 10 Jm , 


1 - exp(-/9t r /i) 


(47) 
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p f _ 

a - o - 180 bars, S = 36 cm and w = 1.6 km. The estimate of G 
is of the same order of magnitude as those obtained by other means 
mentioned earlier. 

A similar technique was employed by Niu (1984/5) to calculate 
the fracture parameters for 49 earthquakes with data compiled by 
Purcaru and Berckhemer (1982). However, he used a slip-weakening 

P ■Jlf 

model which has constant stress a up to 6 = 8 beyond which 
f 

o = o . Such a model does not seem to be in accord with actual 

material behavior. It is easy to show (see, e.g. Rice, 1980) that 

this model reduces the breakdown zone size w in (37) by a factor 

of 4/9. Another method proposed by Ida (1973) and utilized by Aki 

(1979) to obtain fracture parameters for the 1966 Parkfield 

earthquake (Table 4) is also rather similar to the one discussed 

above. The similarity lies in using a source- time model of fault 

P f 

slip to relate the stress-decay (a - a ) to some (indirectly) 

observable time parameter (rise time t described above, and time 

t at which the slip velocity becomes maximum for an in-plane 

shear crack propagating with a uniform velocity used by Aki) and 

to use elastic brittle crack theory and the slip-weakening model 

to estimate G , S and w. While there are some differences in the 
c 

formulae of this class of estimation methods, they are generally 
insignificant when one keeps in mind that the deduced fracture 
parametric values should be regarded as order of magnitude 
estimates only. 

Apart from those described above, Table 4 also summarizes 
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estimates of slip -weakening model and fracture energy parameters 
by Kikuchi and Takeuchi (1970), Husseini et al (1975), Das (1976), 
and a number of other investigators. 

3.4.3 Variations in Fracture Parameters 

It may be seen from the previous discussions and from Table 1 

to 4 that the fracture parameters are quite different between 

those obtained from laboratory tests and those from field 

observations . The laboratory tests for intact rocks give G c on 
4 --2 

the order 10 Jm (and even lower for sawcut rocks, rock joints 

5 8 

and clay) while the field estimates are in the range of 10 to 10 
-2 

Jm (with the exception of Rice and Simons and some of Husseini 's 
estimates) . Similarly the critical slip displacement for 
laboratory samples are in the /^m to mm range, whereas those for 
field estimates are in the cm to m range. These orders of 
magnitude differences are unlikely to be due to temperature or 
normal stress differences since some of Wong's laboratory tests 
were carried out at close to inferred crustal conditions. There 
is reason to believe that natural earthquake faults with 
en-echelon fissures and discontinuities would have "surface 
roughness" orders of magnitude larger than that for laboratory 
specimens, thus contributing to the higher S and G c values. 

To illustrate the plausible dependence of S on surface 
roughness , estimated ranges of 6 from laboratory tests including 
intact and sawcut rocks, jointed rocks and clay, as well as from 
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natural faults are plotted in Fig. 23a. The suggestion is that 

the increasing "surface roughness" of rocks, jointed rocks and 

natural faults are, in addition to normal stress, responsible for 

the increasingly large values of , also shown in Fig. 23b. In 

P f 

summary, these observations indicate that (a - a ) is sensitive 

to normal stress, whereas S is sensitive to surface roughness, and 
- P f 

their product 6 (a - a ) gives the critical energy release rate. 

The laboratory rock data on G for mode II sliding are 

4 -2 

generally of the order 10 Jm . This compares with the much 

lower value in tensile tests which give data mostly in the range 
1 2-2 

of 10 - 10 Jm (see e.g., B. Atkinson, this volume. However, 

1 2 

the 10 - 10 range may be an underestimation of true values; see 

discussion below.) Presumably the micromechanism of the breakdown 
process may be quite different, absorbing much more energy in the 
shear failure mode. 

Lastly, estimation of the size of the breakdown zone w for 

— P f 

laboratory specimens are easily obtained from the S and (a - a ) 
data and applying (37), and they are also listed in Table 1-4. 
These values of w give an indication of minimum characteristic 
dimensions in laboratory specimens for a valid K^ c -test. Wong 
and Rice's analyses of Rummel et al's test results suggest w to be 
in the 0.1- lm range for granitic rocks. Since most laboratory 
specimen sizes are smaller, measurements using standard elastic 
brittle fracture toughness technique are likely to underestimate 
the true K^^-values. For example, the references cited by 
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Atkinson (Table XX in Chapter YY, this volume) give values of 

and Kjj on the order 1 MPa yin. For y - 15-30 GPa, v - 0.25 as 

typical shear modulus and Poisson's ratio for most surface rocks, 

. 2 

this translates to (using (7)) 25-50 Jm . Comparison with 
Table 1 shows that this is at least two orders of magnitude lower 
than that obtained through the slip weakening relation, as 
described in section 3.4.1. 

In general, it would seem advisable to avoid using the 
conventional fracture toughness test, especially for rocks with 
large grain size (and surface roughness, for mode II) unless 
unusually large specimens are used. This statement appears to be 
applicable even for mode' I (tension) fracture toughness testing. 
Figure 24 (after Ingraffea et al, 1984) for example, shows the 
underestimation of fracture toughness for small specimens and 
clearly suggests the size dependence of conventional K^-test 
results. The fitted curves are based on a non-linear analysis to 
be explained in the following section. 


IV. SLIP DISTRIBUTIONS AND INTERACTIONS 

This section describes the representation of slip surfaces 
with generally non-uniform slip in a medium of arbitrary geometry 
and material behavior. In section 4.1, the effect of slippage in 
transferring loads is discussed through Green's functions which 
contain information of the medium material behavior and geometry. 


The formulation results in integral equations when constitutive 
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relations are imposed on the slip surfaces to relate the local 
shear stress to local slippage, as. e.g., in the slip -weakening 
model discussed in Section 3. Such formalism is superior to a 
kinematic description of fault slip because the slip magnitudes 
are part of the solution in solving the problem for a prescribed 
load. The implication is that further physical insight could be 
gained by understanding the slip progression process, which is of 
particular significance in the prediction of slip instabilities. 
Section 4.2 describes the structure of the Green's functions with 
respect to spatial dependence and their homogeneous and 
inhomogeneous parts associated with material boundaries. Selected 
Green's functions most relevant in applications to studies of 
earth faulting are collected in Table 5. For full descriptions of 
such and other Green's functions, the reader should consult the 
references directly. Sections 4.3 and 4.4 review previous studies 
of earth faulting which made use of the methodology described in 
4.1 and 4.2. They are presented in a manner to best illustrate 
the theoretical concepts developed in this and earlier sections . 


4 . 1 Integral Representation and Physical Interpretations 

For a body of any medium with planes of discontinuities , the 
most general representation of the stress state o ^ at a point 
and at time t due to some arbitrary slip S introduced at points Xq 
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and at time t along L (Fig. 25) is 


a . . (x , t) 
ij P 


cr. . (x , t) 
1J P 


- f f- 

J L J-» ^ 


S «(X Q.t ) 

<S p’ S Q' t ' t > d£ dt ' 


(48) 


where a .^.(x^.t) is to be interpreted as the stress state at point 

Xp if no slipping occurs. The stress induced at x^ due to slip at 

all locations x^ is contained in the integral term, carried out 

over all planes of slip displacement discontinuities. The time 

integration is needed in the case of a medium in which memory 

effects are important. These include time -dependent response in a 

viscoelastic medium and diffusion in a fluid- infiltrated 

poro-elastic medium. Indeed the information of the structural 

geometry and the rheology of the stress transmitting medium are 

contained in the Green's function G..(x , x^, t, t ), sometimes 

ij “P “Q 

known as the influence function. It is the fundamental solution 

of the stress a . . at point x at time t due to a unit shear 

~P 

dislocation suddenly introduced at point Xq arid at time t . Many 
elastic and some viscoelastic and poro-elastic solutions for 
various geometries have been obtained. A selection of Green's 
functions useful in describing slippage effects in earthquake 
zones are tabulated in Table 5 . There have been numerous 
interesting applications of these Green’s functions, some of which 
will be described in this section. 
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To give specific interpretation to (48) and to make contact 
with shear rupture models described in Sections 2 and 3, we assume 
here that slip occurs on a single plane (e.g., a boundary between 
two lithospheric plates) between -i<x<i. The plates are subject 
to remote shear loading a° (Fig. 26). a° may be interpreted as 
the shear stress transmitted across a fully locked fault due to 
tectonic scale plate motions. In this case the shear stress 
component in (48) reduces to 


Mx,) 


-f 


i 

-a 


m(i-M 

2n 


V x i 


as(x') 

3x, 


dx. 


(49) 


In (49) the Green's function for an elastic plate under 
generalized plane stress deformation is used (Case 1.2, Table 5). 
This equation has the solution (Muskhelishvili , 1953) 


o° -o(s ) j ds 

(50) 

where "constant" can be determined only from some supplemental 

conditions, e.g. no net dislocation in -i<x^<+i, see (52b) below, 

giving "constant" =0. As a special situation, suppose we impose 

the condition that all points which slip have stress reduced to a 

f f 

constant residual strength a , then a(s) — a in (50) and the slip 


d5( X;L ) 

dx 
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distribution S(x^) may be easily computed in terms of the stress 
drop Aa = <7° - : 


5(x x ) 


. 2Aa /, 2 2 

/»(!+!/) V 1 


(51) 


Equation (51) is, of course, the slip distribution of the shear 
crack model which we described in Section 2 (Eq. (9)) and which we 
used to calculate the energy release rate based on the creep 
displacement rate in central California. The imposition of a 
uniform stress condition on a slip boundary is a characteristic of 
the crack model. The uniform stress condition is appropriate, if 
it is assumed that practically all points on the slip surface have 
undergone sliding exceeding S , in the context of the 
slip-weakening model, or if a quasi -plastic mechanism dominates 
the shear deformation behavior in a narrow shear zone as often 
postulated in the earth’s lower crust where the temperature and 
pressure are high. 

More generally, however, it is appropriate to impose a 
constitutive relation between stress and slip, such as the 
slip-weakening model. In this case, (50) becomes a singular 
integral equation in 8, and 


o(8) 


o 

a 


[ l /*(!+«/) 1 

I , 2 7T X. - X.' 

J -i li 


55(xp 

8x' 


dx. 


(52a) 


the solution of which often requires special numerical methods. 
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Some numerical procedures have been described by Erdogan and Gupta 
(1972), Cleary (1976), Stuart and Mavko (1979) and Fares and Li 
(1986). The equation (52a) and more generally (48) has a unique 
solution only if the net dislocation is specified. For an 
internal discontinuity, a zero net dislocation may be specified as 



as(x') 

<9x,' 


= 0 


(52b) 


It should be noted that (52a) may in fact be regarded as a 
multi-degree-of-freedom system extension of (21) in which the 
stiffness k is analogous to the Green's function here and the 
driving force kS Q is now represented by ct°. 

We now consider the maximum value of a° that can be applied 
to the plate (Fig. 26) before the line of discontinuity -i<x^<i 
extends. This maximum shear load o° may be predicted from brittle 
elastic fracture mechanics if the breakdown zone is much smaller 
than the crack length. Thus from (4a) and (7): 


(a°) - J2p (l+i/)G /n 2 + a ^ 

max c 


(53) 


If the breakdown zone occupies the complete fault then (19) (with 
T identified as ct q ) is applicable, and the maximum allowable load 

p 

would be just the peak value a in the slip -weakening constitutive 
relation. These limits are plotted in normalized form as the 
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slanting and horizontal dashed lines in Fig. 27. The 

characteristic length for normalization of the horizontal 

£ 2 

scale is defined by 2 ^ = 2(l+i/)/iG /(a -a ) which is proportional 
to the breakdown zone size w in (36). For intermediate range of 
the half crack length i, comparable in size to w, (52) together 
with a zero net stress intensity condition (33) was solved 
numerically by Li and Liang (1986) and their result for a linear 
slip-weakening relation is shown as the solid line in Fig. 27. 
(Li and Liang actually solved for the maximum tensile load, but 
the solution coincides with the shear case since the Green's 
functions for the integral equation (52a) are exactly the same for 
both modes of deformation). Clearly, the full numerical solution 
confirms the applicability of the elastic brittle crack model when 
crack size is large in comparison to the breakdown zone size 
(lower left corner of Fig. 27), and the applicability of the 
strength concept at the other extreme (upper right corner of Fig. 
27). 

Conversely, Fig. 27 shows the inadequacy of the strength 
criterion which overestimates the maximum failure load when 
displacement discontinuities exist in the loaded medium. This 
result is consistent with the observation that over-consolidated 
clay slopes often fail by progressive failure at loads much under 
the peak strength of the clay (Bjerrum, 1967). The elastic 
brittle crack model also overestimates the failure load when the 
breakdown zone is comparable in size to the crack length. This 
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observation affords an explanation for the underestimation and 

size -dependence of fracture toughness measured from small 

specimens in the laboratory and based on linear elastic brittle 

crack theory, as described earlier in reference to Fig. 24. 

Indeed the two curve fits in Fig. 24 are based on the predicted 

peak loads from the non-linear analysis shown in Fig. 27 and using 

f 

(4) . (These results hold for both mode I (with a = 0) and mode 
II). To translate from the non-dimensional plot of Fig. 27 to the 

p 

dimensional plot of Fig. 24, we have used a =5 MPa, = 1000, 
1200 psi/in for the lower and upper curves. They are seen to fit 
the experimental data reasonably well. 

4 . 2 Green's Functions and Their Structures 

In the following, we shall further explore the structure of 
the Green's function for dislocation in media with different 
geometries and material properties. The emphasis is to bring out 
the common features between these Green's functions which may on a 
superficial inspection, have little resemblance between each 

t 

other. These common features include the spatial dependence and 
the homogeneous and inhomogeneous parts of the Green's function 
associated with material boundaries. 

For simplicity, we shall confine our focus to 2-D cases, 
although much of the discussions may be extended to 3-D cases as 
well. As further restrictions, we shall limit selected Green's 
functions in Table 5 to static cases, and for geometries most 
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relevant to studying earthquake faulting problems, and then only 
the shear stress component on the line of discontinuity will be 
given. The reader should consult the original literature for a 
full description of the fundamental solutions . We have divided 
the selected Green's functions in Table 5 into three categories: 
they are those for elastic (Case I), viscoelastic (Case II) and 
fluid- inf iltrated poro-elastic media (Case III). Various 

geometries are possible, such as infinite space, half space, plate 
structure, or layered. The dislocation may be of an edge type or 
of a screw type in shear (i.e. in mode II or in mode III). In 
geophysical terminology, the edge dislocation may represent a 
semi-infinite (in length) fault in strike-slip or dip-slip. The 
screw dislocation may represent slip below locked zones in an 
infinite strike-slip fault. The Green's functions in Table 5 are 
given for unit, suddenly introduced dislocations. 

The single major characteristic exhibited by the Green's 
function for all media with different materials and geometries is 
the 1/r singularity (where r is the distance measured from the 
dislocation front). This singular nature of the Green's function 
makes the integral term in (48) a Cauchy principal value integral. 
As mentioned earlier, special numerical techniques are available 
to handle this type of integral. In most cases, terms other than 
the 1/r term occur in the Green's functions. These non-singular 
terms arise because of the presence of material boundaries . For 
example, in a layered medium, these boundaries may divide 
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horizontal regions into layers of different rigidities. Even in a 
half- space, the free -surface exists as a boundary dividing the 
medium into a regular one and one with zero rigidity. The 
non- singular terms are known as the inhomogeneous part of the 
Green's function whereas the singular terms form the homogeneous 
part. In many instances, these inhomogeneous terms have been 
derived by the method of images (see, e.g. Maruyama, 1966, 
Rybicki, 1971) with each term in a summation series representing 
an image point about a plane boundary. As an example, Case 1.4 
shows a screw dislocation at z «■ d below a free surface. The 
homogeneous part of the Green's function is of the form -/i/2?r(z-d) 
with a singularity 'at z - d. The inhomogeneous part is of the 
form p/27r(z+d) from an image source at z = -d, a reflection of the 
primary source about the free surface boundary. 

The elastic rheology assumed for the various geometries shown 
in Case I of Table 4 is plainly an idealization of the mechanical 
behavior of the earth's upper crust. Nevertheless the use of 
elasticity is often justified for the study of short time 
response. The plate geometry (e.g. case 1.1 plane stress) is also 
useful to describe the very long time response of the lithosphere 
when the asthenosphere is fully relaxed, i.e. the lithosphere has 
negligible basal traction and can be therefore treated as a free 
floating plate. 

The viscoelastic behavior of the asthenosphere has been 
suggested to be responsible for many observable time -dependent 
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phenomena. For example, Bott and Dean (1973), Anderson (1975), 

It 

Toksoz et al (1979) , Lehner et al (1981) , and Li and Kisslinger 
(1984/85) have studied the diffusion of stress along a plate 
boundary as a model for migration of great ruptures and for 
filling-in of seismic gaps. Nur and Mavko (1974), Rundle and 
Jackson (1977), Thatcher and Rundle (1979), Thatcher et al (1980), 
Thatcher (1982, 1983, 1984), Melosh and Fleitout (1982), Melosh 

and Raeffky (1983) , Thatcher and Rundle (1984) , Thatcher and 
Fujita (1984), and Li and Rice (1986) have studied the 
time -dependent post-seismic reloading of the lithosphere. Li and 
Rice (1983a, b) have studied the stiffening effect of the 
lithosphere/asthenosphere system as a model for stabilization 
against fault instability and for a precursory period during which 
local plate boundary straining accelerates to failure (see Section 
3.2.2). A common thread of these models is the recognition of the 
coupling between the elastic lithosphere and the viscoelastic 
asthenosphere , with the latter providing the time -dependent 
effect. Usually the time scale of the modelled phenomenon 
provides a rough constraint on the relaxation time and the 
viscosity parameter of the asthenosphere, although other 
non- earthquake related phenomena such as isostatic rebound from 
glacial loading have often been used to estimate the viscosity 
parameter of the asthenosphere. 


Case II. 1 in Table 5 shows a screw dislocation in an elastic 


plate underlain by a Maxwell viscoelastic half-space. Again the 
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Green's function (Bonafede et al , 1986) has a homogeneous and an 
inhomogeneous part due to the presence of boundaries . For the 
short time response, the asthenosphere behaves like an elastic 
body and the stress expression given in Case II. 1 is reduced to 
that in 1.6 with d<H and - ^2 f° r a screw dislocation in an 
elastic half- space. This is the high stiffness limit. (See also 
the discussion in connection with Fig. 14.). For the long time 
response, the asthenosphere is completely relaxed and the stress 
expression given in II. 1. is reduced to that in 1.3 for a screw 
dislocation in an elastic strip, which corresponds to the low 
stiffness limit. As the asthenosphere relaxes between these 
limits, the fault (modelled by the dislocation or a superposition 
of them) is reloaded and the deformation field on the ground 
surface also changes. The time scale for these time -dependent 
transients is given by r *■= 2 rj/p (rj -= viscosity and p = shear 
modulus common to both lithosphere and asthenosphere) in this 
model . 

The above discussion may be extended to case II. 2 which shows 

an edge dislocation suddenly imposed in an elastic plate coupled 

to a viscoelastic foundation through a modified Elsasser model 

(Rice, 1980, Lehner et al, 1981). In the long time limit, the 

system reduces to that of an edge dislocation in a free floating 

plate (Case 1.2, plane stress), as can be shown by taking the 

limit of the time -dependent part of the Green's function, i.e. 

7 (x.. , t-+°°) *= — . For this model, the relaxation time is given by 
j. x n 
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p/a - ( 7 r 2 H/32h)(2r?// i ). 

While there have been several viscoelastic solutions reported 
in the literature, most of them are limited to the time -dependent 
surface displacements only. This is clearly due to the interest 
of modelling post- seismic ground movements based on kinematic 
dislocation models. Unfortunately, Green's functions of the type 
we discuss here are not widely available. 

The time -dependence afforded by a fluid- inf iltrated 

poro-elastic medium comes from the diffusion process associated 

with pore -fluid flow. Such time -dependence has been used as a 

means of modelling after- shock distributions (Booker, 1974, Nur 

and Booker, 1972, Li et al, 1986) and water well level 

fluctuations (Roeloffs and Rudnicki , 1984/85). The stiffening 

effect of an undrained medium on stabilization against faulting 

has been discussed by Rice and Simons (1976) and Rice (1979). 

Cases III.l and 2 show the Green's functions for an edge 

dislocation in such a medium for a permeable (Rice and Cleary, 

1976) and an impermeable (Rudnicki, 1986) fault. The time scale 

2 

for both is controlled by the relaxation time 4c/x^ where c is a 
coefficient of consolidation (see Rice and Cleary, 1976). In the 
long time relaxed limit, the Green's function for Case III.l 
reduces to that of Case 1.1 (plane strain). 

It should be noted that Table 5 represents a small set of 
available Green's functions in the literature. As an example, the 
Green's function for edge dislocation near a circular cavity in an 



75 


elastic medium has been derived by Dundurs and Mura (1964) . Such 
a Green's function could be used to construct the integral 
equation to describe the sliding of a joint near a rock tunnel in 
geotechnical engineering. The possibilities are unlimited, and 
the reader is urged to read the above discussions as a general 
framework which may be specialized to particular applications, 
with the use of the proper Green's function. Many other Green's 
functions can be found in Mura (1982). For many structural 
geometries containing displacement discontinuities, the Boundary 
Element Method coupled with the appropriate Green's function can 
often provide a powerful tool of analysis superior to the Finite 
Element Method (Fares and Li, 1986). 

4. 3 Applications to Dip Slip Faulting 

Two dimensional kinematic models have been used by Chinnery 
and Petrak (1968) and Freund and Barnett (1976) to simulate 
surface vertical movements due to dip -slip faulting on vertical 
and dipping faults. The two dimensionality of these models is 
usually justified on grounds that dip slip faultings occur on 
fault planes with lengths much longer than the widths (in the dip 
direction) . In this section we describe a non-kinematic model due 
to Dmowska (1973) and Dmowska and Kostrov (1973) where the fault 
slip is not preassigned. Instead the fault surface is assumed to 
have uniform and constant shear resistance. This corresponds to 
the assumption of the elastic brittle crack model. Our discussion 
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will begin with the integral representation described in Section 
4.1, and making use of the appropriate Green's function. 

Using the coordinate system shown in Figure 28, and 
specializing to the shear component in the plane of the fault 
surface (48) becomes 


a(s) - a°(s ) - J 2 G(s-s ) 3 ^ S — ds 
S 1 


(54) 


where the time integral has been dropped for the case when the 
material is elastic (with no time -dependent behavior). The 
appropriate Green's function is due to Freund and Barnett (1976) 
(see also Dmowska (1973)) and is shown in Table 5, Case 1.5, 


G(s-s ) 


2m(l 


* 

L s-s 


n \ / * n 5-n 
2 A (s ) s 
n-o n 

2 2 3 

[s +s -2ss cos2a] 


(55) 


where A^ are functions of the dip angle a and are given in Table 5 
(1.5). The presence of the traction free ground surface is 
accounted for by the inhomogeneous part of the Green's function. 
In (54) o° is interpreted as the preexisting tectonic load if the 
fault plane is locked. (In most physical situations, only the 
change in tectonic load is important since only the change, and 
not the true value, in surface deformation can be measured.) 
Assuming that the fault plane slips with uniform shear resistance, 
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the left hand side of (54) is then a constant a , and (54) reduces 
to 


2 G(s-s ) ds (56) 
1 

The singular integral equation (56) has been solved by Dmowska 
(1973) using quadrature formula to discretize the integral, which 
reduces to a set of linear algebraic equations that can be easily 
handled by a computer. The resulting slip distribution, which 
scales with a°-a ^ , is plotted schematically in Figure 28 along the 
fault plane. Using a vertical displacement influence function 
derived by Freund and Barnett, the associated vertical movement on 
the ground surface can also be predicted and is schematically 
sketched in Figure 28 as well. 

The model described above may be made more realistic by 
considering a more sophisticated constitutive relation on the 
fault plane. For example, Dmowska indicated a method of 
incorporating pressure-sensitive friction effects. In that case, 
the shear resistance on the fault would depend on the tectonic 
normal stress acting across the fault and the friction 
coefficient. In addition, slip on the fault plane (in the 
presence of the traction- free boundary) would induce normal stress 
changes and the right hand side of (54) would have an extra term 
similar to the integral term but with a Green's function relating 


^°(s) 


r 
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normal stress change to shear slip. The inclusion of friction 
does not cause any more difficulty in the numerical solution of 
the singular integral equation. 

The analysis could be made even more complete when a 
slip -weakening constitutive law relating shear stress and slip is 
prescribed on the fault plane, as done by Stuart (1979a). The 
si ip -weakening law employed by Stuart (similar to (38)) reflects 
not only a gradual degradation of slip resistance with increased 
amount of slip, but also reflects the increasing ductility (shear 
flow as opposed to brittle fracture) with increasing depth. 
Although Stuart used the finite element method, the formulation 
described above is quite suitable to treat the problem. The only 
difference introduced is on the left hand sides of (54) , where a 
is now made to depend on 6 through the slip -weakening law. 
Introduction of such a term makes the singular integral equation 
non-linear, and the resulting set of non-linear algebraic 
equations will have to be solved by means of an iterative scheme. 

The advantage of such non-kinematic models is that the fault 
slip can come out as part of the solution, and in general is more 
realistic than an imposed uniform dislocation or uniform stress. 
Such non-uniform slip distribution would clearly influence the 
predicted surface deformation behavior, especially if the tip of 
the fault plane is relatively close to the ground surface. In 
addition, it is possible to simulate the progressive failure 
process in response to increased tectonic loading. The failure 
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law is inherent in the slip-weakening model, and Stuart used this 
to calculate the vertical movement at various stages prior to the 
1971 San Fernando earthquake. These results are reproduced in 
Figure 29a and is seen to qualitatively fit the available geodetic 
data. The fault geometry and the variation of peak stress is 
shown in Fig. 29b. It is interesting to note that fault slip 
occurs mostly below 20 km down dip prior to 1969, and a catch up 
process of rapid slip near the hypocenter occurs from 1969 up to 
the 1971 rupture (Fig. 29c). Inferences of such no.n-kinematic 
models provide a means of studying the failure process leading up 
to the slip instability, an earthquake analogue. In general, it 
may be expected that surface deformation may show characteristics 
associated with the approach of an instability, such as 
accelerated vertical movement or strain rates. (See, e.g. Stuart 
and Mavko, 1979 for a detailed discussion of slip instability in 
the context of a strike-slip fault.) Although present available 
data on such precursory signals are scant, precursory deformation 
may nevertheless be useful for assisting future earthquake 
forecasting efforts. 

4 . 4 Application to Slip- stress Interaction Along an 

Inhomogeneous Fault 

Several directly or indirectly observable fault zone 
behaviors suggest that fault surface strength (resistance to slip 
motion) is spatially inhomogeneous. These observations include 
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seismicity concentrations, and episodic creep and repeated 
ruptures of different fault segments. A direct result of fault 
strength inhomogeneity is the non-uniform distribution of fault 
slip and stress accumulation, which no doubt influence surface 
deformation behavior. Thus interpretation of geodetic 

measurements must consider not only depth changes of fault 
property but also along-strike fault property changes, 
particularly where measurements are made close to junctures where 
segments of very different fault behavior occur. In addition, and 
especially relevant to the discussions in this chapter the 
stressing of a locked location must be sensitive to the slippage 
of nearby creeping segments, and such slip-stress interaction must 
be accounted for when considering the processes leading to the 
nucleation of a shear rupture . 

Based on the integral representation described in Section 
4.1, Tse et al (1985) analyzed the stressing of locked patches 
along a creeping fault. Equation (52a), where the Green's 
function for a mode II edge dislocation in an elastic plate (Case 
I. I in Table 5) has been used, describes the mechanics of the 
two-dimensional lithospheric plate shown in Figure 30a. The plate 
is loaded by tectonic stresses a° , and the plate boundary responds 
with a distribution of shear stress ct(x^) and slip £(x^). These 
parameters must necessarily have quantities averaged over the 
thickness of the plate, which is treated as undergoing plane 
stress deformation. It is possible to incorporate the depth-wise 
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change of fault zone properties at the plate boundary by 
considering a cross section at x^, as shown in Figure 30b. This 
section is assumed to undergo antiplane strain deformation, the 
analysis of which provides a spring relationship connecting the 
thickness averaged stress cr(x^) and the thickness averaged slip 
5(x^) through a local spring constant k(x^). It is in this spring 
constant where the details of the depthwise changes in fault 

properties are incorporated. This spring relation is used for the 
left hand side of (52a) which again results in a singular integral 
equation. The formulation of a quasi-three dimensional problem 
described above is really a generalization of the powerful 

line-spring procedure introduced by Rice and Levy (1972) for 
treatment of part- through surface cracks in tension- loaded elastic 
plates or shells. Parks (1981) has shown the remarkable accuracy 
of the approximate procedure in calaulating stress intensity 

factors by comparing it to full scale 3-D finite element 
calculations. The procedure has been used by Li and Rice 

(1983a, b) to analyze strike-slip ruptures in tectonic plates as 
described in section 3.2.2. 

As a specific model, Tse et al (1985) approximated the fault 
zone as sliding under constant stress (taken as zero reference 
stress) at all depths except for a locked seismogenic zone, as 
shown schematically in Figure 4. The free sliding to a depth of 
b(x^) is meant to represent shallow fault creep, and the free 
sliding below the seismogenic zone (between z=b(x^) and z-H-a(x^)) 
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is meant to represent shear flow under essentially constant 
stress. It is useful to note that the geometry of the locked 
patch along- strike is then defined through the dependence of the 
parameters a(x^) and b(x^) on . The assumptions described above 
in essence define an anti-plane strain problem of an elastic strip 
containing two surface edge cracks of depths a and b. The 
solution was obtained by conformal mapping technique and the 
resulting spring constant, which relates the local stress a(x^) to 
the local slip S(x^) in Fig. 30a by a(x^) - k(x^) 5(x^), is (Tse 
et al, 1985) 

r *»(*,) *b(x )-, 

k( V - (fg) / In [2/[cos — + cos— g-ijj (57) 

The stress intensity factors are given in (6) . Note that when a 
and b vanish, the spring constant approaches infinity. This means 
that the local fault slip 5(x^) must be zero for any finitely 
imposed stress when that local segment is fully locked. In (52a) 
the integral limits at -2 and 2 imply that 5=0 beyond this range 
of interest, and this results in a jump in stress a at this 
junction due to the displacement discontinuity. To overcome this 
unnatural artifact, Tse et al modified this assumption to one 
where the stress falls to the tectonic load level a° beyond the 
range -i<x<i, and the assumed uniform slip there can easily be 
computed using (57) for a given far field spring constant k . 
Thus (52a) has an additional term contributed by such far field 
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uniform dislocation, resulting in the form 


kCxpsup + 


r _i_ 3 t(x i ) ) 

J x^ 3Xl dx l 


o fj + H&tKi _l 1 

L . ^. x 2 J 


(58) 


Equation (58) is used to study the stresssing processes of various 
fault locking geometries. Figure 31a shows one of these 
simulations with the lower margin of the locked region chosen from 
depth of the seismicity considerations. For a loading rate <7° of 
0.3 x 10 ^p/yr, the predicted thickness -averaged and surface slip 
rates are shown in Fig. 31b. Geodetic creep data are shown as the 
various symbols in the same figure. In Fig. 31c the thickness 
average stressing rate is shown to vanish in the creep zone and 
falls to the tectonic level at large x^, as required. The 
stressing rate is very high at the tip of the submerged locked 
patch at x^“ 35km, a consequence of interaction between the free 
slip to its left and the sudden locking to its right. Another 
interesting result obtained from the analysis of Tse et al is the 
estimation of the fracture energy release rate of 10 Jm along 
the lower margin of the locked 1857 rupture zone, based on 
required stressing rate to match observed creep data and an 
assumed earthquake cycle time of 150 yr. This order of magnitude 
estimation is again consistent with other estimates already 
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mentioned in section 3 (see also Table 4) . Tse et al noted that 
their attempt to model the locked patch at Parkfield (based on the 
local seismicity and creep data) was limited by the short 
wavelength geometric changes along strike and the basic 
requirement of long wavelength changes in the line -spring 
formulation. 

Stuart et al (1985) solved the same problem using a three 
dimensional version of (48). A three-dimensional Green's function 
due to Chinnery (1963) for a rectangular dislocation patch in an 
elastic half space was employed. Shear resistance on the fault 
plane again takes the form of a bell -shaped slip -weakening law 
(38) . Their model parameters were constrained by repeated 
measurements of fault creep. Figure 32 shows one of their 
computed results compared to creep data near the Parkfield region 
for more than a decade. While an instability event occurring at 
Parkfield is placed at around 1987, Stuart et al cautioned that 
the data would not be sufficient to constrain all the model 
parameters until the fault creep enters the (nonlinear) 
accelerating stage. 

The model of Stuart et al appears tp be more suitable to 
analyzing the Parkfield region because of the inherent 
three-dimensional nature of the patch geometry. For elongated 
patches, such as that recently analyzed by Stuart (1984/5) for the 
500 km segment of the San Andreas fault south of Parkfield, the 
line -spring formulation used by Tse et al should be quite 
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adequate. A characteristic of Stuart and co-worker's models, as 
opposed to many of the available kinematic models in the 
literature, is the ability to analyze the process leading to an 
instability, which might form the corner stone of any earthquake 
forecast model. This is not possible for kinematic models even if 
they 'are able to simulate available surface deformation data, 
since no failure criterion (such as in the form of critical energy 
release rate or slip -weakening law) is employed to track the 
progression of fault slip. 

The line-spring procedure described above in connection with 
Tse et al's work could track the progressive failure process if a 
failure criterion is imposed. For the mode III edge crack model 
(Figure 4) used, the suitable criterion would be a critical energy 
release rate. Li and Fares (1986) studied the stress accumulation 
and slip distribution at the junction of a creep segment and an 
adjacent segment where slip penetration into the seismogenic zone 
occurs under increasing tectonic load. No shallow creep was 
simulated in that study (i.e. b(x)-0 in Figure 4 and in (57)). In 
anticipation of future studies of multiple lines of interacting 
displacement discontinuities, (48) was recast into an indirect 
Boundary Element formulation (Fares and Li, 1986). Following 
Stuart (1979a, b) , but in terms of energy release rate based on 
elastic brittle crack mechanics, the failure criterion was a 
depth- dependent one, as shown in Fig. 17a. An interesting result 
of that analysis was the prediction of a long-term stable slip 
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rate distribution of a parabolic shape in the creep zone, 
decreasing gradually into the adjacent zone which undergoes slip 
penetration and is capable of seismic rupture. With extensive 
slip penetration, the slip-softening behavior becomes more evident 
at this segment, eventually leading to a loss of equilibrium of 
quasi-static slip. This process is accompanied by slip rate 
acceleration, which exceeds that inside the creep zone, as shown 
in Figure 33. 

While the models described in this section and in section 4.3 
incorporate important elements for the study of the instability 
process, results from such instability models nevertheless have to 
be treated with caution. This is because of the assumption of 
pure elastic behavior in the body containing the planes of slip 
discontinuities. In the real earth, the elastic lithosphere is 
underlain by a viscoelastic upper mantle and, possibly, contains a 
viscoelastic lower crust. Time -dependent phenomena attributed to 
the viscoelastic relaxation effect has been described in section 
(4.2). To incorporate the viscoelastic effects, it would be 
necessary to use one of the Green's functions of the type listed 
in Table 5 case II, and the singular integral equation (48) 
requires both a spatial integration and a time integration to 
account for the memory of past slip events. The full solution 
of such an equation is presently not available in earthquake 
instability analyses. A reduced form of (48) where slip is 
averaged over the length of progressive slip zone penetration 
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(e.g, in a seismic gap) has been considered by Li and Rice 
( 1983a, b) (see Section 3.2.2). They found that the instability 
predicted from the models described above corresponds to the onset 
of a quasi-static self-driven period in which stable sliding is 
still maintained but slip acceleration would be inevitable even if 
the tectonic load is kept constant. Physically, this implies that 
a precursory stage, whose time duration is associated with the 

i 

viscosity parameter of the asthenosphere , may precede a dynamic 
instability, or an earthquake rupture. For short term earthquake 
forecasting of great ruptures, it appears to be important to 
capture and interpret the seismic and geodetic data in this 
precursory stage. Future studies of this type of model, with full 
solution of (48), should provide further insight into the time and 
spatial redistribution of surface deformation associated with 
spreading of the softening zone from one or more nucleation points 
on the eventual fault plane prior to a great rupture. 

V. SUMMARY AND CONCLUSION 

This article has focused on the fundamentals of theoretical 
modelling of shear rupture. The slip-weakening model is used as a 
means to unify the discussion, with the elastic brittle crack 
model as one limiting case of extreme non-uniform slip and the 
strength model as the other limiting case of essentially uniform 
slip. The slip -weakening model is regarded as a general 
constitutive law for slip surfaces. Implications in stability of 
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slip systems, and on the extraction of material fracture 
resistance parameters are discussed. Non-kinematic models of 
faulting in various geometries are reviewed and Green's functions 
are described in the context of distributed slip for media of 
elastic, viscoelastic and poro-elastic behaviors. 

Although the available theories of shear rupture have 
provided much insight into understanding the mechanics of earth 
faulting, a complete understanding of many phenomena remains out 
of reach. Recent advances in rate and state dependent 
constitutive laws based on careful experimental observations 
appear to provide a rich foundation on which the transition of one 
earthquake cycle to another could be better understood. There 
appears to be a need to study non-kinematic models in media with 
inelastic behavior in order to understand natural phenomena 
sensitive to the time -dependent rheology of the earth. Finally, 
natural faults are never ideally straight and' with only a single 
strand, and fault surfaces are likely to have mechanical 
properties varying along- strike and with depth. These 
characteristics call for 3-D modelling, in order to describe 
non-uniform distributed slip on multiple non-linear fault strands. 
It is likely that the fault constitutive laws, fault geometry and 
the medium rheology all play important roles in controlling the 
time and location of the nucleation of a slip instability. 
Advances in understanding slip rupture behavior will have to come 
from laboratory and in- situ experimentation’ and from analytic and 
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numerical modelling. While progress has been and will continue to 
be made in the near future in these directions, it may be expected 
that serious obstacles exist. For example, it is not clear how 
one might translate experimental observations in the laboratory to 
the field given the orders of magnitude difference in the values 
of slip -weakening model parameters obtained in the laboratory and 
those estimated from field observations. Also, the more 
sophisticated and complete a model is, even if sufficient 
computational power is available (which is not necessarily the 
case for non-linear problems of the type suggested by the rate and 
state dependent friction laws), the more model parameters will 
need to • be constrained. At the present time, available 
geophysical data, especially those collected precursory to a large 
plate boundary rupture, is extremely limited. 

VI . ACKNOWLEDGEMENTS 

The author acknowledges useful discussions with the following 
individuals: R. Dmowska, H. Einstein, N. Fares, H.S. Lim, J.R. 
Rice and W.D. Stuart. H.S. Lim contributed significantly to the 
preparation of the tables and figures in this text. I thank 
C. Benoit for typing this manuscript under a very tight schedule. 
This article was written when my child Dustin was born. I am 
grateful to my wife, Stella, for her immaculate patience, to my 
mother-in-law, for her helpful assistance at home, and to my son 
for sleeping through the night. This work would not have been 



90 


completed without the generous cooperation of all the individuals 
mentioned here. 

This article was prepared under support of the NSF Geophysics 
Program, the USGS Earthquake Hazards Reduction Program, and the 
NASA Crustal Dynamics Program. 


91 


VII. REFERENCES 

Aki, K. (1979). J. Geophys . Res. 84, pp. 6140-6148. 

Aki , K. and Richards, P.G. (1980). Quantitative Seismology Theory 
and Methods, Vol I & II, W.H. Freeman, San Francisco. 

Allen, C.R. and Smith, S.W. (1966), Bull Seism. Soc . Amer. 56, 
pp. 966-967. 

Anderson, D.L. (1975). Science, 187, pp . 1077-1079. 

Ashby, M. and Hallam, S.D. (1986). Acta MetaI1.34, 3, pp. 497-510. 

Barenblatt, G.I. (1962). In "Advances in Applied Mechanics", 
Academic Press Inc., London, V.7, pp. 55-125. 

Barton, N. (1972). Int. J. Rock Mech. Min. Sci., 9, pp. 579-602. 

Bjerrura, L. (1967). Trans. Am. Soc. Civil Engrs . SM 93, pp. 3-49. 

Bonafede, M. , Boschi, E. and Dragoni , M. (1986). J. Geophys. 

Res . , in press . 

Booker, J.R. (1974). J. Geophys. Res. 79, pp. 2037-2044. 

Bott, M.H.P., and Dean, D.S. (1973). Nature, 243, pp . 339-341. 

Brune , J.N. (1970). J. Geophys. Res. 75, pp. 4997-5009. 

Burford, R.O. and Harsh, P.W. (1980). Bull. Seism. Soc. Am. 70, 
pp. 1233-1261. 

Cherepanov G.P. (1968). Int. J. Solids Struct. 4, pp. 811-831. 

Chinnery, M.A. (1961). Bull. Seismol. Soc. Amer. 51, pp. 355-372. 

Chinnery, M.A. (1963). Bull. Seismol. Soc. Amer. 53, pp. 921-932. 



92 


Chinnery, M.A. (1970). In "Earthquake Displacement Fields and the 
Rotation of the Earth" (edited by L. Mansinha et al) . Dordrecht 
Reidel, pp. 308. 

Chinnery, M.A. and J.A. Petrak (1967). Tectonophysics 5, 
pp. 513-529. 

Cleary, M.P. (1976), International J. for Numerical Methods in 
Engineering , 10, pp . 679-720. 

Coulson, J.H. (1972). In Proceedings of XIII 1971 U.S. Symposium 
on Rock Mechanics, (E.J. Cording, ed.), pp. 77-105. 

Das, S. (1976). In "A Numerical Study of Rupture Propagation and 
Earthquake Source Mechanism", PhD thesis, M.I.T.. 

Dieterich, J.H. (1978). Pure Applied Geophys . 116, pp. 790-806. 

Dieterich, J.H. (1979). J. Geophys. Res. 84, pp. 2161-2168. 

Dmowska, R. (1973). Publ. Inst. Geophys. Acad. Sci. 62, 
pp. 124-139. 

Dmowska, R. and Kostrov, B.V.,(1973). Archives of Mechanics 25, 3, 
pp . 421-440, Warszawa. 

Dmowska, R. and Rice, J.R. (1986). In "Continuum Theories in Solid 
Earth Physics" (R. Teisseyre, ed.) pp . 187-255, Elsevier 
Publishing . 

Dugdale, D.S. (1960). J. of Mech. Phys . Solids 8, pp. 100-104. 

Dundurs and Mura (1964). J. Mech. Phys. Solids 12, pp. 177-189. 

Erdogan, F. and Gupta, G.D. (1972). Quarterly Applied Mathematics , 


pp. 525-534. 



93 


Eshelby, J.D. (1957). P roc. Roy. Soc. (London) Ser. A 241, 
pp. 376-396. 

Evans, J.B. and Wong, T.F. (1985). In "Mechanics of Geomaterials, 
Rock Concrete and Soil", (Z. Bazant ed.).pp. 189-210. John Wiley 
& Sons, Chichester. 

Fares, N. and Li. V.C. (1986), J. Engineer. Frac Mech., in press. 

Freund, L.B. (1976). Mechanics Today 3, Nemat-Nasser (ed.) S., 
Pergamon, N.Y., pp . 55-91. 

Freund, L.B. and Barnett, D.M. (1976a). Bull. Seism. Soc. Amer . 

66, 3, pp . 667-675. 

Freund L.B. and Barnett, D.M. (1976b). Bull. Seism. Soc. Amer. 66, 
pp. 2083-2084. 

Goodman, R.E. (1970). In "Determination of the In Situ Modulus of 
Deformation of Rock", ASTM STP 477, pp. 174-196. 

Griffith, A. A. (1920). Phil. Trans. Roy. Soc. London, Ser. A221, 
pp. 163-198. 

Head, A.K. (1953). Proc . Phys . Soc. (London) B66, pp . 793-801. 

Hillerborg, A., Modeer, M. and Petersson P.E. (1976). Cem. and 
Conor. Res. 6, pp . 773-782. 

Hirth, J.P. and Lothe, J.L. (1986). Theory of Dislocations , McGraw 
Hill, NY. 

Husseini, M.I., Jovanovich, D.B, Randall, M.J. and Freund, L.B. 
(1975). Geophys. J.R. Astr. Soc. 43, pp. 367-385. 

Ida, Y. (1972). J Geophys Res. 77., pp . 3796-3805. 

Ida, Y. (1973). Bull. Seismol. Soc. Amer. 83 , pp . 959-968. 



94 


Ingraffea, A.R. , Gunsallus, K.L. , Beech, J.F. and Nelson, P.P. 

(1984). ASTM STP 885, pp. 152-166. 

Irwin, G.R. (1960). In "Structural Mechanics: Proc. of 1st Sym. 
on Nav. Struc. Mech, 1958" (Goodier and Hoff, eds) pp. 

557-591 , Pergamon Press, NY. 

•Jaeger, J.C. and Cook, N.G.W. (1969). "Fundamentals of Rock 
Mechanics", Methuen and Co., Ltd., London. 

Kanamori, H. and Anderson, D.L. (1975). Bull Seism Soc . Amer . 

65, pp. 1073-1095. 

Kasahara, K. (1981). "Earthquake Mechanics", Cambridge University 
Press, Cambridge. 

Kikuchi, M. and Takeuchi, H. (1970). Zisin, 23, pp. 304-312. 

King, C.Y. , Nason, R.D. and Tocher, D. (1973). Phil. Trans. Roy. 

Soc. London, Ser. A, 274, pp. 255-360. 

Kostrov B.V. and Nikitin L.V. (1970). Arch. Mech. Stos., 22, 750. 
Lehner, F.K., Li, V.C. and Rice, J.R. (1981). J. Geophys . Res. 

86, pp. 6155-6169. 

Li, V.C. and Fares, N. (1986). In "Earthquake Source Mechanics" 
5th Maurice Ewing Symposium, AGU, in press. 

Li, V.C. and Kisslinger, C. (1984/85). Pure and Applied Geophys. 
122, pp. 812-830. 

Li, V.C. and Liang, E. (1986). J. of Engineering Mechanics, ASCE, 
in press. 

Li, V.C. and Rice, J.R. (1983a). J. Geophys. Res. 88, pp . 


4231-4246. 



95 


Li, V.C. and Rice, J.R. (1983b) Bull Seismo. Soc . Amer . 73 
pp. 1415-1434. 

Li, V.C. and Rice, J.R. (1986). In preparation. 

Li, V.C., Seale, S.H. and Cao, T. (1986) Tectonophysics . 

Lisowski, M. and Prescott, W.H. (1981). Bull. Seism. Soc. Am. 71, 
pp. 1607-1624. 

Maruyama, T. (1963). Bull. Earthq. Res. Inst., Tokyo U. 41, 
pp . 46 - 86 . 

Maruyama, T. (1964). Bull. Earthq. Res. Inst., Tokyo U. 42, 
pp. 289-368. 

Maruyama, T. (1966). Bull. Earthq. Res. Inst., Tokyo U., 44 
pp. 811-71. 

Melosh, H.J. and Fleitout, L. (1982). Geophys . Res. Lett. 9, pp . 
21-24. 

Melosh, H.J. and Raeffky, A. (1983). J. Geophys. Res. 88, pp. 
515-526. 

Minster J. B.and Jordan, T.H. (1984) . in "Tectonics and 

Sedimentation along the California Margin" (Crouch, J.K. and 
Bachman, S.B., eds), Pacific Section S.E.P.M., vol. 38, pp . 

1-16. 

Morgenstern, N.R. and Tchalenko, J.S. (1967). Geotechnique 17 
pp. 309-328. 

Mura, T. (1982). "Micromechanics of Defects in Solids", Martinus 


Nijhoff, London. 



96 


Muskhelishvili , N.I. (1953) "Singular Integral Equations", 
translated by J.R.M. Radok, Noordhoff, Groningen. 

Nabarro, F.R. N. (1952). In "Advances in Physics, Vol. 1, 
pp. 169-179. 

Nemat-Nasser , S. and Horii, H. (1982). J. Geophys Res. 87, 
pp. 6805-6821. 

Niu, Z. (1984/5). Pure and Applied Geophys. 122, N.5, pp. 645-661. 

Nur, A. and Booker, J.R. (1972). Science, 183, pp . 204-206. 

Nur, A. and Mavko, G. (1974). Science, 183, pp. 204-206. 

Okubo, P. and Dieterich, J.H. (1984). J. Geophys. Res. 89, 
pp. 5817-5827. 

Palmer, A.C. and Rice, J.R. (1973). Proc . Roy. Soc . Lond . , Ser. A 
332, 527-548. 

Papageorgiou, A.S. andAki, K. (1983). Bull Seis. Soc. Am. 73, 
pp. 953-978. 

Parker, A.P. (1981). "The Mechanics of Fracture and Fatigue", 

E. & F.N. Spon Ltd. , London. 

Parks, D.M. (1981). Trans. ASME, J. Pressure Vessel Technol. 103, 
pp. 246-254. 

Paris, P.C. and Sih, G.C. (1965). In "Fracture Toughness Testing 
and its Applications", STP 381, pp. 30-76, ASTM , Phila. 

Purcaru, G. and Berckhemer, H. (1982). Tectonophysics , 84, 
pp. 57-128. 

Rice, J.R. (1968a). In "Fracture, An Advanced Treatise", 

Vol. 2, H. Liebowitz (ed.), Academic Press, NY, pp 191-331. 



97 


Rice, J.R. 

(1968b) 

. J. Appl. Mech. 35, pp . 379 

-386. 

Rice , 

, J.R. 

(1979) . 

Gerlands Beitr. Geophysik, 

Leipzig 88, 2 

S. 

91-127 

• 



Rice , 

, J.R. 

(1980) . 

In "Physics of the Earth's 

Interior" 


(A.M. Dziewonski and E. Boschi, eds) , pp. 555-649, Italian 
Physical Soc., No. Holland. 

Rice, J.R. (1983). Pure and Applied Geophys . 121. No. 3, 
pp. 443-475. 

Rice, J.R. (1984). Proc. of 1st Inter. Congress on Rockbursts 
and Seismicity in Hines, Johannesburg, 1982, SAIMM, (Gay, N.C. 
and Wainwright, eds.), pp. 57-62. 

Rice, J.R. and Cleary, M.P. (1976). Rev. Geophys. Space Phys . , 24, 
(2), pp. 227-241. 

Rice, J.R. and Levy, N. (1972). J. Appl . Hech. 39, 185-194, 1972. 

Rice, J.R. and Simons, D.A. (1976). J. Geophys. Res. 81, 

5322-5334. 

Roeloffs, E. and Rudnicki, J.W. (1984/85). Pure and Appl. Geophys. 
122, 2-4, pp. 560-582. 

Rudnicki, J.W. (1980). Ann. Rev. Earth Planet Sci. 8, pp. 489-525. 

Rudnicki, J.W. (1986). In "Earthquake Source Mechanics", 5th 
Maurice Ewing Sym. , AGU, in press. 

Ruina, A. (1983). J. Geophys. Res. 88, pp . 10359-10370. 

Rummel, F. , Alheid, H.J. and Frohn, C. (1978). Pure and Applied 
Geophys. 116, pp. 743-764. 

Rundle, J.B, (1986) J Geophys. Res., Vol. 91, pp. 1947-1959. 



98 


Rundle, J.B. and Jackson, D.D. (1977), Geophys . J.R. Astro. Soc . 
49, pp. 575-591. 

Rybicki, K. (1971). Bull Seis. Soc. Amer . , 61, pp. 79-92. 

Rybicki, K. (1986). In "Continuum Theories in Solid Earth Physics" 
(R. Teisseyre, ed.) pp . 18-186, Elsevier, NY. 

Sacks, I.S., Suyehino, S. Linde, A.T. and Snoke , J.A. (1978). 
Nature 275, pp. 599-602. 

Sano, 0. (1981). In "Fundamental Study on the Mechanisms of 

Brittle Fracture of Rock", D. Eng. Thesis, Kyoto Univ. , Kyoto, 
Japan. 

Savage, J.C. and L.M. Hastie (1966). J. Geophys. Res. 71, 
pp .• 4897-4904. 

Schulz, S.S., Mavko, G.M., Burford, R.O. and Stuart, W.D. (1972). 

J. Geophys. Res. 87, pp . 6977-6982. 

Sieh, K.E. (1984). J. Geophys. Res. 89, pp . 7641-7670. 

Skempton, A.W. (1964). Geotechnique 14, pp. 77-101. 

Spottiswoode , S.M. (1980). J. Geophys. Res., submitted. 

Stekekee, J.A. (1958a). Can J. Rhys., 36, pp. 192-205. 

Stekekee, J.A. (1958b). Can J. Phys . , 36, pp. 1168-98. 

Stuart, W.D. (1979a). Science, 203, pp. 907-910. 

Stuart, W.D. (1979b). J. Geophys. Res. 84, pp . 1063-1070 . 

Stuart, W.D. (1984/5). Pure and Applied Geophys. 122, pp. 793-811. 
Stuart, W.D., Archuleta, R.J. and Lindh, A.G. (1985). J. Geophys. 
Res. 90, pp . 592-604. 


(7 - SL 



99 


Stuart, W.D. and Mavko, G.M. (1979). J. Geophys . Res. 84 
pp. 2153-2160. 

Tada, H. , Paris, P.C. and Irwin, G. (1973). In "The Stress 

Analysis of Cracks Handbook", Del Research Corp., Hellertown, 
Pa. 

Thatcher, W. (1982). Nature 299, pp . 12-13. 

Thatcher, W. (1983). J. Geophys. Res. 88, pp. 5893-5902. 
Thatcher, W. (1984). J. Geophys. Res. 89, pp. 1867-1873. 
Thatcher, W. and Fujita, N. (1984). J. Geophys. Res. 89, pp. 
3102-3106. 

Thatcher, W. , Matsuda, T., Kato, T. and Rundle, J. B. (1980). 

J. Geophys. Res. 85, pp . 6429-6439. 

Thatcher, W. and Rundle, G.B. (1979). J. Geophys. Res. 84, pp. 
5540-5556. 

Thatcher, W. and Rundle, G.B. (1984). J. Geophys. Res. 89, pp . 
7631-7640. 

II 

Toksoz, M.N. , Shakal, A.F. and Michael, A.J. (1979). Pure Appl. 
Geophys., 117, pp . 1258-1270. 

Tse, S.T., Dmowska, R. and Rice, J.R. (1985). Bull. Seism. Soc . 
Arne r 75, 3, pp. 709-736. 

Tse, S.T. and Rice, J.R. (1986). J. Geophys. Res., in press. 
Tullis, T. and Weeks, J. (1986). In "Earthquake Source 
Mechanics", 5th Maurice Ewing Symposium, AGU, in press. 
Turcotte, D. and Spence, D.A. (1974). J. Geophys. Res. 79, 


pp. 4407 -4412. 



100 


Volterra, V. (1907). Ann. Sci. Ecole Norm. Super., Paris 24, 
pp. 401-517. 

Walsh, J.B. (1969). J. Geophys. Res. 74 (8), pp. 2070-2080. 
Wawersik, W.R. and Brace, W.F. (1971). Rock Mech . , pp. 61-85. 
Wong, T.F. (1982a). J. Geophys Res. 87, pp. 990-1000. 

Wong, T.F. (1982b). Mechanics of Materials 1, pp. 3-17. 

Wong, T.F. (1986). In "Earthquake Source Mechanics", 5th Maurice 
Ewing Symposium, AGU, in press. 

Yip, C.K. (1979). In "Shear Strength and Deformability of Rock 
Joints", M.I.T. MS Thesis 



101 


Table 

Table 

Table 

Table 

Table 

Table 


TABLE CAPTIONS 


la: Slip-weakening Model Parameters for Intact Rocks, 
lb: Slip -weakening Model Parameters for Sawcut Rocks. 

2 : Slip -weakening Model Parameters for Rock Joints (all at 

room temperature) . 

3 : Slip -weakening Model Parameters for Over-consolidated 

Clay . 

4 : Slip-weakening Model Parameters for Natural Faults. 

5 : Green's Functions for Dislocations in Elastic, 

Viscoelastic and Fluid- inf iltrated Poro-elastic 


Media. 
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FIGURE CAPTIONS 


Fig. 1 


Fig. 2 


Fig. 3 


Fig. 4 


Fig. 5 


Tip of a shear crack in mode II and mode III deformation. 
The darker shade denotes an annular region in which the 
asymptotic crack tip stress fields given by (1) and (2) 
are valid. 

Center crack with half crack length in an elastic plate 

loaded by remote shear stress a° . Crack faces have 

f 

uniform shear resistance a . 

Semi- inf inte crack in mode II deformation in an elastic 
infinite body. Crack faces are loaded by line forces P 
(per unit thickness) at a distance b from crack tip. 

Double edge cracks with crack lengths a and b in mode III 
deformation in an elastic strip loaded at the remote 
boundaries by thickness averaged stress a. The thickness 
averaged slip displacement S is that associated with the 
presence of the cracks. 

Schematic illustration of the extension process of the 
crack tip and the associated work absorbed by relaxing the 
stress (la) with simultaneous crack tip displacement 


(lb). 


Fig. 6 Slip rate data from a 200 km trace of the San Andreas 
fault in central California (after Burford and Harsh, 
1980; Lisowski and Prescott, 1981; Schulz et al , 1972). 
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The curve fit is from the elastic brittle center-crack 
model (9) , with coordinates origin set at 10 km west of 
Monarch Peak. The broad- scale geodolite data carries 
stronger weight because of the 2-D (thickness -averaged) 
nature of (9) . 

Fig. 7 Tube of material cut out by contour T near crack tip, to 
illustrate the balance of energy flux into this tube of 
material and energy absorbed by elastic work in A, 
frictional work on L, and energy which drives crack 
extension, G. 

Fig. 8 The J- integral applied to a mode II shear crack, to relate 
the energy release rate G to J and the frictional 
dissipation. 

Fig. 9 Application of the J- integral to the San Andreas fault to 
extract the critical energy release rate G c associated 
with the 1857 Ft. Tejon earthquake M - 8.3. 

Fig. 10 (a) Triaxial test results for initially intact granite 
samples reported in Rummel et al (1978) for three 
different confining pressures, (b) Slip-weakening 
branches deduced by Rice (1980, 1984) from raw data in 
(a). 

Fig. 11 Slip-weakening curves for four types of rock joints (from 
Goodman, 1970). Note that when slip -hardening occurs as 
in case 4b, the slip -weakening model and associated 
theories are not applicable. 
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Fig. 12 Schematic plot of the constitutive slip-weakening model, 

with (a) showing the post-peak weakening with S, and the 

rigid unloading branches, and (b) showing the increase of 

peak strength and residual frictional strength as a 

function of effective normal stress, and (c) showing the 

decrease of peak strength and residual frictional 

strength as a function of temperature. Note that (b) and 

(c) have been drawn to illustrate the general increase 

P f 

followed by decrease of strength drop (ct -a ) with a 
and the general decrease in strength drop with T. 

Fig. 13 (a) A single -degree -of- freedom spring-block model, with 
loading through imposed displacement 5 q , and load 
transmitted through a spring with stiffness k. (b) Trace 
of equilibrium loads a , a , .... a , and corresponding 

A d hi 

slips 5^, ...,5g. Instability sets in at E, when 

equilibrium cannot be maintained, followed by slip 
acceleration and rapid stress drop rate approaching 
infinity, as illustrated in (c) . Reestablishment of 
equilibrium can be at any of the points F, G or H. (d) 

For a stiffer spring and the same slip-weakening 
relationship, the unloading lines are steeper and no 
dynamic instability occurs. The stress may drop and the 
slip may accelerate as in (e) , but their time rate of 
change do not approach infinity. 
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Fig. 14 (a) A single-degree-of-freedom system loaded through a 

standard viscoelastic element, (b) The point I represents 
initial instability, when the loading system is fully 
relaxed, with stiffness k^. The point D represents 
dynamic instability when the loading system reaches its 
maximum stiffness . Between I and D, the block is 

self-driven . 

Fig. 15 (a) Upward progression of slip in a seismic gap zone in 
the elastic lithosphere of thickness H underlain by a 
viscoelastic foundation of thickness h. The lithosphere 
is treated as undergoing 2-D plane stress deformation as 
shown in (b) which depicts loading of the plate by a° , 
and with (thickness averaged) stress a and slip 5 in the 
seismic gap zone. The o-S relation at the plate boundary 
is derived from the anti -plane strain mode III elastic 
brittle crack model shown in (c) . The resulting 2-D 
problem is further reduced to a single-degree-of-freedom 
model by assuming that a is uniform in the seismic gap 
zone . 

Fig. 16 Time -dependent compliance of the coupled-plate system 

(solid line), and its approximation by a single parameter 
standard linear model (dashed line) for two seismic gap 
lengths . 

Fig. 17 (a) Assumed fracture energy variation with depth. 

(b) Thickness -averaged stress versus slip {o-S) relation 
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for elastic-brittle crack model of the antiplane rupture 
progression shown in Fig. 16c. Based on (18) and the 
fracture energy distribution as in (a). The point P 
represents the peak stress state and the point M 
represents the maximum slip possible for the o-S 
relation shown. Instability must occur between P and M. 

Fig. 18 Solution of (24) for (a) crack penetration, (b) averaged 
stress, and (c) averaged slip, as a function of time t. 
The subscripts P, I, and D for the normalizing parameters 
denote the peak, initial instability, and dynamic 
instability states. Each plot shows the solution 
corresponding to two loading rates R. 

Fig. 19 (a) Stress and slip distributions near a crack tip with a 
breakdown zone in which the deformation behavior is 
governed by the slip -weakening relation. The weakening 
branch is shown in (b). 

Fig. 20 (a) Assumed linear variation of stress in breakdown zone 
(34) and (b) the corresponding si ip -weakening relation, 
for estimation of the breakdown zone size w. After 
Palmer and Rice (1973). 

Fig. 21 (a) Rock specimen loaded triaxially. (b) Experimental 
output of differential stress versus axial shortening. 

(c) Relation between axial relative movement of sliding 
surfaces and slip S . (d) Derived slip-weakening curves 
from (b) . From Rice (1980). 



107 


Fig. 22 Composite of critical energy release rate versus normal 
stress, from two test series on San Marcos gabro and 1 
Fichtelbirge granite. After Wong (1986). 

Fig. 23 (a) Nominal slip displacement ranges and (b) critical 
energy release rates for various geo-materials, from 
experimental testing and inferences from field 
observations . 

Fig. 24 Size dependence of apparent mode I fracture toughness of 
various rock types on crack length in laboratory 
specimens. Data from Ingraffea et al (1984). The curve 
fits are based on numerical solution of the non-linear 
singular integral equation (52) by Li and Liang (1986), 
assuming a linear slip-weakening relation. The upper fit 
is for a plateau value (large crack length a) of of 
1200 psi J in, and the lower fit is for a plateau value of 
Kj c of 1000 psi Jin. See text for further details and 
also refer to Fig. 27. 

Fig. 25 A general body containing a line of displacement 

discontinuity L. Slip S at x induces stress a . . at x . 

q ij “P 

Fig. 26 Elastic plate containing a single line of displacement 
discontnuity with stress distribution ct(x^) and slip 
distribution 6(x^). The plate is loaded remotely by o° . 

Fig. 27 Prediction of maximum applied load for the center crack 
plate structure shown as insert, from elastic brittle 
crack theory (slanted dashed line) and from strength 
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criteria (horizontal dashed line) . The solid line is 
predicted by including a breakdown zone where material 
deforms according to a linear slip -weakening relation, 
and is numerically obtained by solving (52) . 

Fig. 28 Shear slip distribution and induced vertical ground 
displacement in dip slip faulting from s^ to s^. 

Fig. 29 (a) Comparison of observed and predicted uplift at various 
times prior and up to the 1971 San Fernando earthquake, 
based on a slip-weakening fault model with geometry shown 
in (b) . The fault slips for various time periods are 
shown in (c) . 

Fig. 30 (a) Elastic model of tectonic plate assumed .to undergo 

plane stress deformation, and loaded by remote stress o ° , 
with stress and slip distributions a and 6 at plate 
boundary, (b) A local section of the plate boundary, 
assumed to undergo anti-plane strain deformation. The 
shaded fault zone can be modelled by any appropriate 
constitutive law, and the resulting relation between 
ct(x^) and 5(x^) defines the local spring constant k(x^) . 

Fig. 31 (a) One of several geometries of locked patches on a fault 
surface analyzed by Tse et al (1985) for stressing of the 
Parkfield region of the San Andreas fault. The lower 
margin of the patch has been chosen from depth of the 
seismicity consideration, (b) Comparison of model 
prediction based on the solution of (58) to geodetic slip 
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rate data, (c) Computed thickness averaged stress 
distribution along strike. From Tse et al (1985). 

Fig. 32 (a) Map view of the San Andreas fault strand near 

Parkfield. (b) Geometry of locked patches, (c) Comparison 
of model prediction of fault creep at various locations 
shown in (a) with creepmeter data. The theoretical creep 
has been multipled by 0.8 to compensate for 
underestimation of fault slip measured by creepmeters. 

From Stuart et al (1985). 

Fig. 33 Computed slip rate as a function of distance along strike 
at various stages prior to instability, in a creep zone 
centered at x^*=+2H. Outside the creep zone, upward 
penetration of the brittle zone occurs with increasing 
tectonic loading. Note that the parabolic slip rate 
distribution qualitatively agrees with contemporary data, 
as shown in Fig. 31b, at the three earlier time steps 
shown. As instability approaches, the slip rate at the 
edge of the creep zone accelerates and exceeds that inside 
the creep zone. From Li and Fares (1986) 
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Table lb : Slip-weakening Model Parameters for Sawcut Rocks 
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Table 2 : Slip-weakening llodel Parameters for Pocl; Joints (all at room temp.) 
















Table 3 : Slip-weakening Model Parameters for over-consolidated clay 
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Table 4 : Slip-weakening Model Parameters for Natural Faults 
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CASE . STRUCTURAL GEOMETRY AND RHEOLOGY HOMOGENEOUS INHOMOGENEOUS 
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CASE STRUCTURAL GEOMETRY AND RHEOLOGY HOMOGENEOUS REFERENCE 
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• ALINEMENT ARRAY (66-79) 
■ CREEPMETER (68-80) 

A SHORT-RANGE (75-79) 

-f GEODOLITE (69-78) 



DISTANCE ALONG TRACE (km) 






Figure 







DISPLACEMENT, cm DISPLACEMENT, cm 



0 1 2 3.4 


DISPLACEMENT, cm 



type type of surfaces typically exhibiting behavior in this class 

1 HEALED joints and incipient joints 

2 CLEAN, SMOOTH FRACTURES 

2a POLISHED 

2b UNPOLISHED (ROUGH SAW CUT) 

3 CLEAN, ROUGH FRACTURES 

3a ARTIFICIAL EXTENSION FRACTURES AND DISTURBED SAMPLES 
3b UNDISTURBED SAMPLES 

4 FILLED JOINTS, SHEARED ZONES, SHALE PARTINGS, AND SMOOTH BEDDING 

4a DRY OR SLIGHTLY MOIST 
4b WET; THIN SEAM 
4c WET; THICK SEAM 


Figure 11 








Figure 13 












AVERAGE STRESS AT PLATE BOUNDARY 
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Figure 20 







SAWCUT INTACT 0/P ROCK NATURAL 

ROCK ROCK CLA.Y JOINTS FAULTS 
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Figure 26 
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Figure 33 



